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The  activation  of  spoilers  on  the  upper  surface  of  a  wing  is  a  relatively  new  method  of 
achieving  longitudinal  and  lateral  control  of  ram-air  canopies,  often  referre 
oarafoih  No  numerical  studies,  however,  have  fully  investigated  the  3-D  aerodynamic 
performance of 'these  bleed-air  actuators.  Simulation  results  are  presented  for  a  finite 
span  ram-air  canopy  geometry  and  several  configurations  amenable  for  comparison 
wdth  wind  tunnel  experiments  and  flight  test  data.  Assessments  arc  also  made  between 
the  resulting  flow-fields  and  previously  reported  2-D  computational  fluid  j  wmiis 
fCFD)  results.  Evaluations  include  cases  with  an  asymmetrical  trailing  e  Se  L 
to  simidate  a  turn  and  the  inclusion  of  two  different  bleed-air  vents.  Lift,  drag  and  roll 
coefficients  correla  ie  well  w  ith  wind  tunnel  and  flight  test  results. 


Nomenclature 


AR 

b 

aspect  ratio,  b/S 
span 

L/D 

Re 

c 

airfoil/wing  chord 

t 

cL 

lift  coefficient,  L/qrx  S 

cD 

drag  coefficient.  D  j  qccS 

s 

c, 

rolling  moment,  L/ q^Sb 

V 

Cm 

pitching  moment,  M  / qa,Se 

u * 

+ 

C„ 

yawing  moment,  N  j qyS h 

y 

cr 

side  coefficient,  Y /  qrrS 

a 

D 

drag  force 

P 

L 

lift  force 

V 

lift  to  drag  ratio 
Reynolds  Number,  Vc/v 
Time 

non-dimensional  time,  t*Vfc 

Area 

Velocity 

friction  velocity  at  wail 
y-plus,  u*y/v 
angle  of  attack 
Density 
air  viscosity 


L  Introduction 


An  autonomously  guided  cargo  deliver  consists  of  a  ram-air  canopy,  suspension  lines  connecting  he 
canopy  t“bome  guidance  unit  which  in  turn  is  connected  to  a  payload.  The  original  ram-a,r  inflated 
wine  was  proposed  by  Domina  Jalbert  in  1961 '  with  the  trailing  edge  being  deflected  down  for  Uteral  an 
longitudinal  flight  control.  An  asymmetric  deflection  provides  turning  authority  skid  steering,  such  that  a 
rfoht  turn  is  produced  by  the  right  side  deflection.  A  symmetric  deflection,  on  the  other  hand  changes  the 
descent  rat/and  glide* angle.  For  both  the  asymmetric  and  symmetric  tail  deflection  the  resulting 
dynamics  has  led  to  the  term  “brakes”  being  the  common  vernacular  to  identify  this  control  mechanism. 
Ware  and  Hassel”  conducted  wind  tunnel  tests  of  various  ram  air  canopies  and  note  symmetric  deflection 


*  Senior  Research  Aerospace  Engineer,  A1AA  Senior  Member 
t  Senior  Aerospace  Engineer,  AIAA  Senior  Member 
t  Research  Associate,  AIAA  Senior  Member 
5  Research  Fellow,  AIAA  Senior  Member 
'  Director,  AIAA  Senior  Member 
I'  Professor  of  Aeronautics,  AIAA  Associate  Fellow 
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the  unrestrained  case,  and  6°,  “when  maximum  flap  was  applied  ”  ’  f°r 
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parachutists.  '  "  feh  performance  improvements  to  be  adopted  by  personnel 
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NASA10  anH  dl-n^11  ^an<?Pje^  t0  autonomously deliver  payloads  began  in  earnest  in  the  1990s  bv 
range  of  10  lbs  to  42Klbs12  Both  th^  v  .  TDAnc  &C  parachute  systems  across  a  payload  weight 

Higgins  documents  the  potential  use  of  a  bleed  air  vent  system  for  autonomously  o  7  i  ■  . 
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comprehensive  Verification  and  Validation  (V&V)  initiated  by  Natick  Soldier,  Research,  Development, 
and  Engineering  Center  (NSRDEC)20'21.  In  addition,  this  study  extends  2-D  bleed-air  vent  simulations  to 
3-D.  Results  are  presented  using  the  finite  volume  method  solvers  Kestrel  3  and  Cobalt 4,  which  have 
been  extensively  validated  for  unsteady  CFD  for  fixed  wing  aircraft.  Kestrel  is  a  fixed  wing  multi¬ 
physics  computational  engineering  product  of  the  High  Performance  Computing  and  Modernization 
Program’s  (HPCMP)  Computational  Research  and  Engineering  for  Acquisition  Tools  and 
Environments-Air  Vehicle  (CREATE 1  '"-A V)  program.  The  results  also  support  parallel  work  in 
Ghoreyshi  et  al25  which  specifically,  addresses  meshing  considerations  for  these  types  of  coupled  cavity- 
airfoil  geometries.  The  paper  first  reviews  the  grid  generation,  flow  solver,  and  computational  parameters. 
Test  cases  are  then  described,  and  results  are  presented.  Finally,  conclusions  and  recommendations  are 
given. 


II.  Computational  Setup 

As  previously  mentioned,  this  computational  effort  was  conducted  in  parallel  to  the  results  reported  in 
Ghoreyshi  et  al25.  Details  of  the  grid  convergence  study  for  various  geometries  as  well  as  validation 
results  are  included  in  that  presentation. 

A.  Geometry 

The  baseline  geometries  used  for  the  validation  simulations  are  shown  in  Figure  1.  These  geometries 
are  derived  from  the  cut  pattern  used  in  the  manufacturing  process  of  a  ram-air  canopy  being  drop  tested 
and  which  is  designed  to  enable  precision  landings  with  a  nearly  vertical  final  descent  profile.  The  solid 
models  represent  a  half  span  of  the  actual  wing,  such  that  the  model  aspect  ratio,  AR,  is  1 .  I  hey  were 
created  in  a  modular  manner  to  enable  the  construction  of  multiple  wind  funnel  configurations  that  would 
be  used  for  validation  and  testing.  Four  attachment  holes  for  the  wind  tunnel  mounting  can  be  seen  in  the 
root  of  each  wing.  To  provide  a  well-posed  airfoil  shape  for  the  rounded  closed  inlet  case,  a  cubic  spline 
was  used.  The  rounded  leading  edge  may  be  removed,  as  emphasized  by  the  red  circles  in  Figure  1 ,  to 
produce  a  closed,  or  open,  flat  nose  with  an  angle  the  same  as  the  open  inlet  ram-air  canopy.  Though  the 
rounded  and  fiat  nose  leading  edge  geometries  do  not  represent  a  ram-air  parachute,  they  are  critical 
models  for  conducting  simulation  validation,  and  they  serve  as  a  reference  for  understanding  the  flow 
characteristics  of  the  ram-air  canopy’s  hybrid  cavity-airfoil  design.  An  additional  configuration  of 
interest  included  the  50%  asymmetric  deflection  of  the  trailing  edge  of  the  airfoil,  highlighted  by  the 
green  circles  in  Figurel.  Finally,  Figure  Id.  illustrates  the  combination  of  the  leading  and  trailing  edge 
modifications. 

Future  rigid  geometries  and  solid  models  will  include  additional  features  that  more  accurately  capture 
steady  flight  configurations  of  the  ram-air  canopy.  These  characteristics  include  the  inflated  deformation 
of  the  ram-air  canopy  cells  as  well  as  the  arc  anhedral.  Open  inlet  geometries,  and  meshes,  with  a  straight 
and  50%  deflected  trailing  edge  were  computationally  modeled  and  simulations  completed,  but  no  solid 
models  were  constructed  for  these  specific  configurations.  The  open  inlet  simulation  data  will,  however, 
be  compared  against  future  wind  tunnel  experiments. 

Eslambolchi  and  Johari26  conducted  rigid  canopy  simulations  of  a  different  baseline  airfoil,  with 
inflation  and  arc  anhedral  included,  to  assess  the  fidelity  of  the  lifting-line  theory  in  predicting  the  forces 
on  the  canopy  during  steady  glide.  The  authors  found  the  forces  on  the  canopy  during  steady  glide  using 
finite  wing  theory  to  be  only  accurate  to  approximately  30%  for  the  expected  operational  range  of  angle 
of  attack  for  their  canopy.  These  differences  were  attributed  to  a  variety  of  factors  which  accentuate  the 
difficulties  of  high-fidelity  ram-air  modeling  and  simulation  as  well  as  obtaining  accurate  aerodynamic 
force  and  moment  measurements  from  parachute  drop  tests. 

Using  the  current  geometries,  two  additional  configurations  including  two  different  versions  of  the 
bleed-air  vents  were  simulated.  The  placement  of  the  vents  was  restricted  to  the  middle  of  the  ram-air 
airfoil,  and  followed  the  designs  for  which  previous  2-D  CFD  simulations  had  been  completed.21'  Vent 
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width  was  approximately  14%  of  the  half  span.  This  size  was  chosen  to  model  the  equivalent  vent  size  on 
an  airdrop  canopy.  Finally,  an  open  inlet  geometry  with  an  asymmetric  TD  of  50%  was  simulated,  and 
represent,  to  the  authors  knowledge,  the  first  such  simulations.  These  results  will  be  presented  and 
compared  against  2-D  CFD  flowfields,  and  as  mentioned  earlier,  they  will  also  be  compared  against  wind 
tunnel  experiments  in  a  future  effort. 

B.  Computational  Grid 


Both  structured  and  hybrid  grids  tor  a  half  wing  with  a  symmetry  plane  were  tested,  and  will  be 
identified  in  the  presentation  of  results.  The  "structured"  meshes  (structured  everywhere  except  the  wing 
tip  surface  and  the  volume  mesh  from  the  tip  to  the  far-field  boundary)  had  between  15M  and  39M  cells, 
and  were  created  using  Pointwise.  The  hybrid  meshes  used  a  refined  boundary  layer  grid  along  the  airfoil 
surfaces  to  capture  the  boundary  layer  while  a  tetrahedral  grid  was  used  for  the  far  field.  The 
unstructured  meshes  ranged  from  5.5M  to  5.9M  cells,  and  they  were  created  using  CREATE™-AV 
Capstone.  Sample  projections  of  the  various  meshes  used  in  these  simulations  are  illustrated  in  Figure  2. 

Except  for  the  BA  grids,  the  grids  have  a  y+  <  1  which  ensures  the  first  computational  cell  above  the 
wall  is  located  within  the  viscous  sub-layer.  Figure  3  illustrates  the  wall  y+  for  the  open  inlet-no  TD 
configuration.  For  the  BA  vent  grids,  they*  around  the  vent  flaps  is  approximately  2.  The  growth  rate  in 
the  viscous  layer  is  <  1 .25  for  all  grids.  For  the  symmetry  wall  the  grid  spacing  was  0. 1 .  Additionally,  the 
farfield  was  located  at  +/-  10c. 

C.  Flow  Regime  and  Simulation  Parameters 

CFD  results  were  computed  using  both  the  commercially  available  code  Cobalt  and  the  DoD- 
developed  solver  within  the  Kestrel  software  suite.  Both  codes  are  based  on  finite-volume  methods  and 
in  addition  to  performing  traditional  Reynolds  Averaged  Navier-Stokes  (RANS)  calculations,  also  have 
the  capability  to  perform  Delayed-Detached  Eddy  Simulations  (DDES).  Cobalt  simulations  were  run  on 
the  Cray  XE6  at  Engineering  Research  Development  Center  (ERDC)  DoD  Supercomputing  Resource 
Centers  (DSRC)  while  Kestrel  was  used  on  the  SGI  ICE  X  System  located  at  the  AFRL  DSRCr 

The  numerical  method  used  by  Cobalt  is  based  on  Godunov's  first-order  accurate,  cell-centered,  finite 
volume,  exact  Riemann  solution  method  applicable  to  arbitrary  cell  topologies27.  The  spatial  operator  uses 
the  exact  Riemann  Solver  of  Gottlieb  and  Groth,  least  squares  gradient  calculations  using  QR 
factorization  to  provide  second  order  accuracy  in  space,  and  TVD  flux  limiters  to  limit  extremes  at  cell 
faces  ‘ .  A  point  implicit  method  using  analytic  first-order  inviscid  and  viscous  Jacobians  is  used  for 
advancement  of  the  discretized  system.  For  time-accurate  computations,  a  second  order  accurate  Newton 
sub-iteration  scheme  is  employed.  Parallel  performance  is  achieved  using  the  ParMETIS  domain 
decomposition  library  for  optimal  load  balancing  with  a  minimal  surface  interface  between  zones39.  The 
code  uses  Message  Passing  Interface  (MPI)  for  communication  between  processors,  with  parallel 
efficiencies  above  95%  on  as  many  as  1024  processors  and  scales  linearly  up  to  4000  processors. 

The  flow  solver  component  of  Kestrel  solves  the  unsteady,  three-dimensional,  compressible  Navier- 
Stokes  equations  on  hybrid  unstructured  grids.  Its  foundation  is  also  based  on  Godunov's  first-order 
accurate,  exact  Rjemann  solver.  Second-order  spatial  accuracy  is  obtained  through  a  least  squares 
reconstruction.  I  he  code  uses  an  implicit  Newton  sub-iteration  method  to  improve  time  accuracy  as  well. 
Kestrel  receives  an  extensible  Markup  Language  (XML)  input  file  generated  by  Kestrel  User  Interface 
and  stores  the  solution  convergence  and  volume  results  in  a  common  data  structure  for  later  use  by  the 
Output  Manager  component.  Some  available  turbulence  models  are  the  Spalart-Allmaras  (SA)  model,  SA 
for  Rotation/Curvature  (SARC),  and  DDES  with  SARC.  For  these  simulations  the  inviscid  flux 
calculations  were  computed  with  the  Harten,  Lax,  van  Leer,  and  Einfeldt  (HLLE)  Riemann  solver,  and 
the  Barth/Jespersen  flux  differencing  scheme  was  used. 

Table  1  documents  the  parameters  and  conditions  used  for  the  unsteady,  2nd  order  in  time,  simulations 
which  resulted  in  a  Reynolds  number.  Re  —  1 .4  x  10’’  with  a  non-dimensional  time  step  A /*=  0.007.  Given 
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these  values,  Figure  4  shows  typical  convergence  histories  for  the  lift,  drag,  pitch,  roll,  and  yaw 
coefficients.  These  specific  plots  are  for  the  FB  (flat  leading  edge,  50%  TD)  geometry.  Convergence 
occurs  between  1000  and  1500  iterations. 


Table  1.  General  Job  File  Inputs  and  “Test”  Matrix 


Parameter  or  Property 

Value  or  Condition 

Gas  Model 

ideal  Gas 

Turbulence  Model 

SARC-DDES 

Newton  Subiterations 

3 

Iterations 

5000 

Mach 

0.25 

a  (deg) 

0, 4,8, 12, 16,20 

Static  Pressure  (Pa) 

81,850 

Static  Temperature  (K) 

290 

Time  step  (s) 

0.000025 

Velocity  (m/s) 

84 

IV.  Simulation  Results  and  Discussion 

In  this  section,  computational  results  are  presented  for  the  rigid  configurations.  The  first  data  serve  as 
preliminary  validation  of  computational  techniques  and  processes.  Though  simulations  were  conducted 
for  0°  <  a  <  20°,  analyses  of  the  flowfields  will  not  be  discussed  for  all  conditions.  The  set  of  presented 
results  will  illustrate  flowfield  evolution  for  open  and  closed  inlet  geometries  at  a  =  8°  as  this  value  is 
representative  of  the  nominal  angle  of  incidence,  and  thus  angle  of  attack,  used  for  ram-air  systems  used 
for  personnel  and  cargo  airdrop.  Finally,  data  will  be  discussed  for  various  open  inlet  designs  including 
the  straight  trailing  edge,  deflected  trailing  edge,  and  front-  and  rear-flap  bleed-air  vents. 

A.  Validation  with  Wind  Tunnel  Experiments,  RS  Geometry 

In  an  earlier  validation  effort  Ghoreyshi  et  ai.21  compared  CFD  results  against  data  from  experiments 
performed  at  Kyushu  Institute  of  Technology  of  Japan30  with  inconclusive  results.  The  Japanese 
experiment  measured  a  setup  with  a  flexible  wing  surface  supported  by  fixed  ribs  at  a  Mach  number  of 
0.04  and  Reynolds  number  of  2  *  105.  However,  the  simulations  assumed  a  rigid  structure  throughout, 
and  it  is  believed  this  structural  difference  is  the  main  contributor  to  the  discrepancies  among  the 
measured  forces  and  moments.  NSRDEC  and  the  USAFA  subsequently  designed  a  new  set  of 
experiments  using  the  geometries  described  above,  and  a  comprehensive  report  of  the  experimental 
results  will  be  reported  later. 

The  Ci  and  CD  from  the  Kestrel  simulations  and  the  airfoil  wind  tunnel  experimental  values  are 
depicted  in  Figure  5  and  show  good  agreements.  It  should  also  be  noted  that  a  grid  convergence  study, 
which  included  using  Cobalt  was  conducted  as  part  of  these  analyses.  As  a  result  a  grid  size  of 
approximately  6M  cells  was  used  in  Kestrel  for  the  remaining  computations.  Figure  6,  illustrates  the 
attached  flow  about  a  finite  RS  geometry,  and  shows  the  simulation  capturing  the  strong  wing  tip  vortex. 
With  these  validation  results,  the  remainder  of  the  presentation  will  address  the  aerodynamic 
characteristics  of  the  different  geometries  at  a  -  8°. 

B.  Closed  Geometries:  RS,  FS,  RB,  and  FB 

In  Figure  7  a  plot  of  the  isosurface  of  normalized  Q  =  0.1,  i.e.,  the  second  invariant  of  the  gradient 
tensor31,  colored  by  pressure  is  used  to  track  vortical  structures  in  the  flow.  The  predominant  effect  due 
to  geometry  changes  is  seen  in  corresponding  changes  in  the  topology  of  the  wing  tip  vortex.  Three 
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vortex  structures  are  identified  in  Figures  7a.  and  7b, —  tip  nose  ,  tip  upper  surface,  and  tip  lower  surface 
vortices.  These  structures  may  also  be  identified  by  the  reader  in  Figures  7c.  and  7d. 

The  tip  nose  vortex  of  the  RS  geometry  evolves  from  the  separation  of  the  flow  along  a  relatively 
large  portion  of  the  leading  edge  at  the  wing  tip.  The  locally  unsteady  flow  organizes  into  a  structure  at 
.1c  and  follows  the  pressure  gradient  to  the  top  of  the  wing  near  ,25c.  At  that  position  the  tip  nose  vortex 
passes  over  the  tip  upper  surface  vortex  which  has  been  created  near  ,25c.  Additionally,  a  similar  lower 
surface  tip  vortex  has  been  created  at  approximately  .  1 25c,  and  slowly  grows  spatially  while  also 
following  the  pressure  gradient  towards  the  upper  surface  until  it  crosses  the  tip  at  .75c.  For  the  FS  wing, 
the  pointed  nose  at  the  wing  tip  causes  a  concentrated  vortex  to  form  but  with  a  u'eaker  amount  of  energy 
as  compared  to  the  RS  tip  nose  vortex.  As  such,  the  majority  of  the  FS  tip  nose  vortex  joins  with  the 
upper  surface  tip  vortex.  The  remaining  wing  tip  separated  flow,  below  the  pointed  nose,  bursts  into 
several  small  vortical  structures  until  the  majority  of  the  enstrophy  organizes  itself  with  the  lower  surface 
vortex  and  follows  the  favorable  pressure  gradient  around  the  wing  tip  to  .75c  as  did  the  RS  case.  A  small 
vortical  structure  is  seen  to  join  the  upper  surface  vortex  at  ,5c.  The  resulting  lift  is  reduced  and  drag 
increased  for  the  FS  relative  to  the  RS  configuration,  and  the  L/D  is  reduced  by  approximately  4%. 

Figures  7c.  and  7d.  illustrate  the  effect  of  the  asymmetric  trailing  edge  deflection  on  the  closed 
geometries.  In  particular,  the  magnitudes  of  the  lift,  drag,  pitch,  roll,  and  yaw  coefficients  increase,  but 
where  the  RS  and  FS  L/D  are  approximately  7.1  and  6.8  respectively,  the  L/D  for  RB,  5.6,  and  for  BF, 
5.5,  are  significantly  reduced.  Indeed,  the  TD  leads  to  significant  changes  in  all  aerodynamic 
characteristics  with  a  60%  increase  in  CD  (L/D  decreases  by  20%)  and  a  100%  increase  in  C„.  These 
characteristics  align  with  the  use  of  the  trailing  edge  deflection  for  longitudinal  and  lateral  control  of  the 
ram-air  canopies — glide  slope  changes  result  from  symmetric  TD  and  a  turn  is  initiated  by  an  asymmetric 
trailing  edge  deflection  on  the  either  side  of  the  canopy.  The  drag  and  yaw  prevail  over  lift  and  roll 
increases.  Airdrop  tests  of  different  ram-air  geometries  have  shown  control  inputs  with  lower  amounts  of 
TD  which  result  in  initially  opposite  lateral  and  longitudinal  responses,  but  these  conditions  were  not 
investigated  in  this  work.  Table  2  documents  the  changes  in  the  force  and  moment  coefficients  for  the 
FS,  RB,  and  FB  configurations  relative  to  RS. 

Table  2.  Percent  Change  Forces  and  Moments,  Closed  Geometries 

(Relative  to  RS*) 


Coefficient 

FS  (RS) 

RB  (RS) 

FB  (RS) 

CL 

-0.8% 

27.8% 

26.2% 

cD 

3.4% 

63.5% 

65.0% 

L/D 

-4.1% 

-21.8% 

-23.5% 

C, 

4.8% 

40.6% 

45.6% 

cm 

0.9% 

14.3% 

14.7% 

Ci 

-0.6% 

39.0% 

37.1% 

C„ 

8.4% 

95.4% 

101.9% 

A  flow  description  similar  to  the  non-deflected  trailing  edge  cases  follows  for  the  TD  cases,  though  in 
Figure  7d.  at  ,5c  additional  bursting  along  the  outward  face  of  the  wing  tip  has  now  joined  with  the  upper 
surface  vortex.  The  more  favorable  pressure  gradient  for  the  deflected  trailing  edge  cases  also  leads  to  a 
widening  in  the  vortical  structures  at  the  trailing  edge.  Additionally,  the  increased  curvature  at  the  trailing 
edge  induces  separation* 

C.  Open  Inlet  Geometries:  Open,  50%  TD,  and  Bleed-air  Vents 

At  a  near  8°,  except  for  the  rear-flap  bleed-air  vent  configuration,  the  2-D  simulations  of  ram-air 
geometries 1 0  exhibited  a  strong  separation  bubble,  reference  Figures  8  and  9.  In  the  case  of  the  rear-flap 
bleed-air  vent,  the  flow  along  the  upper  surface  of  the  geometry  was  significantly  accelerated  resulting  in 
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an  attached  flow  and  an  order  of  magnitude  increase  in  the  ram-air  airfoil's  L/D.  However,  the  results 
presented  below  support  the  need  to  complete  3-D  simulations  in  order  to  obtain  valid  data  tor  these 
otherwise  geometrically  simple  cavity-airfoil  flows. 

1,  0%and50%TD 

The  percent  change  relative  to  the  corresponding  closed  geometry  for  the  0%  and  50%  TD 
configurations  are  documented  in  Table  3.  Overall  aerodynamic  performance  is  diminished.  In  both 
cases  the  L/D  is  reduced  by  the  same  percentage,  as  compared  to  the  closed  rounded  geometries.  The 
major  differences  in  coefficients  are  in  the  side  force  and  yaw.  For  the  0%  TD  case,  the  percent  relative 
change  for  both  C,.  and  C„  is  twice  as  large  as  the  50%  TD  relative  change,  and  indicates  the  dominant 
contribution  of  the  TD  as  opposed  to  the  change  from  the  closed  to  open  inlet.  In  Figure  10,  the  upper 
surface  flow  is  attached  except  for  the  small  separation  region  near  the  root  of  the  wing.  The  flowfield 
has  a  three-dimensional  structure  similar  to  the  closed  inlet  geometries  RS  and  FS,  but  with  a  reduced  L/D 
=  6.2  and  an  increase  of  over  25%  in  C„. 

Table  3.  Percent  Change  Forces  and  Moments,  Open  Geometries 


(Relative  to  equivalent  closed  geometries.) 


Coefficient 

0%  TD  (RS) 

50%  TD  (RB) 

CL 

-6.3% 

-8.0% 

cD 

8.4% 

6.4% 

L/D 

-13.6% 

-13.5% 

Cv 

10.7% 

5.4% 

Cm 

-6.0% 

-9.3% 

Cl 

-6.6% 

-7.9% 

c„ 

26.2% 

12.1% 

Figure  1 1  shows  a  much  stronger  separated  flow  at  the  wing  root. 

2.  0%  TD,  Front  Flap  and  Rear  Flap  Bleed-air  Vents 

In  Figures  10-13,  the  ram-air  airfoil  surfaces  have  been  illustrated  with  a  transparency  effect  such  that 
0  structures  may  be  seen  in  the  interior  of  the  cells.  While  the  magnitudes  of  the  vortical  structures  are 
significantly  less  than  those  seen  on  the  outer  surface,  there  is  a  nonzero  variation  throughout.  Actual 
ram-air  canopies  typically  include  crossports  between  adjacent  cells,  positioned  at  several  ehordwise 
locations,  to  allow  pressure  equalization  across  the  span  of  the  ram-air  airfoil.  The  crossport  features  are 
important  to  include  for  accurately  modeling  future  fluid  structure  interactions  (FSI),  but  they  are  not 
modeled  in  these  results. 

Table  4.  Percent  Change  Forces  and  Moments,  Bleed- Air  Vents 


(Relative  to  the  0%TD  geometry.) 


Coefficient 

Front  Flap  (0%TD) 

Rear  Flap  (0%TD) 

CL 

-14.1% 

-0.6% 

Cd 

49.8% 

0.6% 

L/D 

-42.6% 

-1.1% 

Cr 

-9.5% 

-1.9% 

Cm 

-21.9% 

-0.8% 

Q 

-12.3% 

-0.4% 

c„ 

59.9% 

1.0% 
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Qualitatively,  the  3-D  BA  vent  flows  mirror  the  behavior  of  the  2-D  simulations  in  Figure  9.  However 
for  the  front  flap  BA,  even  though  the  intersection  of  a  streamwise  plane  centered  at  the  mid-span  point  of 
the  vent  with  the  flow  would  show  a  field  similar  to  the  2-D  field  of  Figure  9a.,  the  3-D  structure  (ref 
Figure  12)  has  complex  transverse  recirculation  components.  The  strength  of  the  front  flap  BA  jet  causes 
a  separation  bubble  which  affects  75%  of  the  ram-air  airfoil  span.  Alternatively,  the  flow  field  for  the 
rear  vent  BA  shows  a  relatively  mild  perturbation  as  compared  with  the  laminar  flow  of  the  no-vent  ram- 
air  geometry.  The  3-D  effect  is  much  less  than  that  predicted  by  2-D  simulations,  and  the  result  also 
follows  more  closely  with  field  data. 

1  able  4  quantitatively  shows  the  relative  effects  of  the  BA  vents  as  compared  to  the  no-vent,  ram-air 
geometry.  For  the  front-flap  BA  vent,  the  decrease  in  lift  and  large  increases,  with  respect  to  the  0%  TD 
baseline,  in  drag  and  yaw  reflects  the  effectiveness  of  the  "spoiler.'’  The  3-D  rear  flap  BA  vent  forces  and 
moments,  though,  indicate  a  significantly  different  overall  behavior  as  compared  to  2-D. 

3.  0%  TD,  50%  TD,  and  Front  Flap  Bleed-air  Vent 


The  data  in  fable  5  identifies  the  different  aerodynamic  contributions  to  longitudinal  and  turning 
response  changes  depending  on  the  type  of  actuation  used.  No  attempt  has  been  made  to  identify 
conditions  for  equivalent  turning  response  or  longitudinal  rate  changes  between  the  two  control  methods. 
While  turning  response  comparisons  require  a  more  detailed  analysis,  for  the  defined  geometries,  the  data 
indicates  the  front  flap  BA  vent  is  a  more  effective  mechanism  for  longitudinal  control  than  TD. 

Table  5.  Percent  Change  Forces  and  Moments 


(Relative  to  the  0%TD  geometry*) 


Coefficient 

50%  TD  (0%  TD) 

Front  Flap  BA  (0%TD) 

CL 

25.5% 

-14.1% 

Co 

60.4% 

49.8% 

L/D 

-21.7% 

-42.6% 

Cy 

33.8% 

-9*5% 

c,„ 

10.2% 

-21.9% 

c, 

37.1% 

-12.3% 

a 

73.6% 

59.9% 

The  positioning  of  the  BA  vent  was  somewhat  arbitrary,  though  the  current  design  lends  itself  to 
modular  application  across  the  surfaces  of  the  ram-air  geometry.  Therefore,  a  geometric  superposition  of 
vents  on  the  surfaces  of  the  canopy  may  be  easily  accommodated.  Trailing  edge  deflection  modifications 
are  less  systematically  employed.  Additionally,  changes  in  the  baseline  airfoil  geometry  will  have 
significant  influence  on  responses  of  each  control  method. 

V.  Conclusion 

Computational  results  for  3-D  ram-air  geometries  including  trailing  edge  deflection  and  bleed-air 
vents  have  produced  a  rich  data  set  which  may  be  used  for  control  law  design.  The  validation  of  the 
numerical  models  against  wind  tunnel  data  suggests  meshes  on  the  order  of  6M  cells  are  sufficient  for 
engineering  analysis  lor  these  conditions,  the  results  confirm  the  application  of  symmetric  trailing  edge 
deflection  for  longitudinal  control  using  a  downward  deflection  on  the  order  of  45‘.  However,  a  single 
BA  vent  produced  twice  the  change  in  L/D  as  was  produced  by  the  TD.  Skid  steering  was  also  verified  as 
the  dominant  turning  response  for  ram-air  canopies,  while  BA  vents  were  shown  to  be  effective  control 
methods  in  this  regard.  Future  work  will  continue  the  validation  of  the  computational  methods  with  wind 
tunnel  experiments  and  airdrop  tests.  New  efforts  will  also  address  BA  vent  location  and  TD/BA  vent 
coupling. 
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a)  Round  Leading  Edge-Straight 
Trailing  Edge  (RS) 


b)  Flat  Leading  Edge-Straight 
Trailing  edge(FS) 


c)  Round  Leading  Edge- Asymmetric  Trailing 
Edge  Deflection  (RB) 


d)  Flat  Leading  Edge-  Asymmetric  Trailing 
Edge  Deflection  (FB) 


Figure  1:  Solid,  Closed  Inlet  Ram-air  Canopy  Geometries 


Figure  2:  Example  Meshes 


V  Pb,5 


Figure  3  :  Mesh  Quality  Example — Open  Inlet  No  TD  (OS) 
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Turne  rul  Tt«e|fcl  Tjrttejif 

Figure  4;  Example  convergence  histories — FB  at  a  =  8°. 


Figure  5:  CFD-to-Wind  Tunnel  comparison  for  the  RS  model 


Figure  6:  Attached  flow  about  a  finite  ram-air  closed  inlet  wing 
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Upper  Vortex 
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Figure  7:  a  =  8°  Simulation 


a,  0%  TD  b.  50%  TD 

Figure  819:  2-D  CFD  Ram-air  Simulations,  a  =  8-5 


13  of  14 


American  Institute  of  Aeronautics  and  Astronautics 


a.  Front-flap  Deflection 


b.  Rear-flap  Deflection 


Figure  91’:  2-D  CFD  Bleed-air  Simulations,  a  =  8.5 


Figu  re  10:  0%  TD,  L/D  =  6.2  Figure  1 1 :  50%  Asymmetric  TD,  L/D  =  4.8 
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Figure  12:  BA  Front-flap  Deflection,  L/D  =  3.5  Figure  13:  BA  Rear-flap  Deflection,  L/D  =  6.1 
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This  paper  provides  an  overview  on  the  grid  quality  and  resolution  effects  on  the  aero¬ 
dynamic  modeling  of  parachute  canopies.  The  CFD  simulations  were  performed  using  the 
Cobalt  flow  solver  on  two  dimensional  canopy  sections  with  open  and  closed  inlets.  Previ¬ 
ous  simulation  results  of  these  geometries  showed  that  the  grid  independence  is  achieved 
for  the  closed  and  open  airfoils  with  grids  containing  around  half  a  million  and  two  million 
cells,  respectively.  Previous  grids  were  either  hybrid  with  prismatic  layers  near  the  wails 
or  multi-block  structures  using  algebraic  grid  generators.  The  results  presented  in  this 
work  show  that  the  grid  independence  of  both  geometries  can  be  achieved  with  even  much 
coarser  grids.  These  grids,  however,  were  generated  with  good  smoothness,  wall  orthog¬ 
onality,  and  skewness.  Cobalt  solver  reports  an  overall  and  a  local  quality  value  for  each 
grid  cell  The  results  show  that  this  quality  value  is  mainly  related  to  the  grid  smooth¬ 
ness  but  does  not  depend  on  the  mesh  skewness  or  the  wall  orthogonality.  Although  a 
smooth  grid  improves  the  Cobalt  quality  value  and  therefore  the  solution  convergence,  but 
it  does  not  always  lead  to  an  accurate  solo f ion.  For  example,  the  unstructured  meshes  with 
anisotropic  cells  near  the  wall  have  very  good  grid  quality  in  Cobalt,  but  they  have  the 
worst  accuracy  between  the  considered  meshes  because  of  the  mesh  poor  skewness  at  the 
walls.  The  results  also  showed  that  in  comparison  to  the  closed  inlets,  the  open  geometry 
solutions  are  less  sensitive  to  the  initial  grid  spacing  and  number  of  constant  spacing  layers 
at  outside  wails.  Also,  the  open  inlet  solutions  do  not  change  with  inside  mesh  resolution 
and  type. 


Nomenclature 


lift  coefficient,  L/q^S 
drag  coefficient,  D/q^Sc 
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D 

drag  force,  N 

GQ 

grid  quality 

GR 

growth  rate  in  the  viscous  layer 

L 

lift  force,  N 

l 

far  field  length  away  from  the  airfoil,  m 

M 

Mach  number,  V/a 

N, 

the  number  of  layers  of  constant  spacing  normal  to  all  viscous  walls 

goc 

dynamic  pressure,  Pa,  pV2  / 2 

to* 

kth  grid  quality  metric  at  ith  cell 

Re 

Reynolds  number,  pVc/p 

t 

time,  s 

V 

freestream  velocity,  m/s 

X,  y ,  z 

grid  coordinates 

Greek 

AsJ 

the  initial  spacing  at  all  viscous  walls 

Q 

angle  of  attack,  rad 

P 

density,  kg/ m3 

P 

air  viscosity 

I.  Introduction 

The  accuracy  and  expediency  of  computational  fluid  dynamic  (CFD)  solutions  depend  not  only  on  the 
underlying  numerical  methods,  but  also  on  the  grid  generation  process,  in  which  the  computational  domain  is 
discretized  into  distinct  sub-domains.  These  sub-domains  are  called  cells  or  grid  blocks.  The  numerical  error 
lises  if  the  physical  geometry  is  not  exactly  represented  by  the  grid.1  Also,  truncation  error  and  machine 
round-off  error  are  always  present  in  the  numerical  solutions  of  partial  or  ordinary  differential  equations,2 
These  errors  depend  on  the  grid  size  such  that  for  a  finer  mesh,  the  truncation  errors  decreases,  but  the 
rounding  errors  will  increase.  Finally,  inaccurate  interpolation  of  the  discrete  solutions  between  grid  blocks 
will  increase  the  numerical  error. 

The  uniform  grids  have  many  advantages  over  non-uniform  grids  including  a  better  accuracy  and  faster 
convergence.  However,  to  resolve  boundary  layer  and  wake  regions  in  viscous  flows,  a  uniform  grid  ap¬ 
proach  would  require  a  very  large  number  of  mesh  points,  which  increases  the  computational  cost.  Practical 
turnaround  times  lor  obtaining  a  CFD  solution  and  availability  of  computer  facilities  would  limit  the  uni¬ 
form  grid  applicability  for  three-dimensional  (3D)  problems.  The  total  number  of  mesh  points,  and  therefore 
computational  cost,  can  be  reduced  by  using  a  non-uniform  grid  approach,  in  which  the  grid  points  are 
clustered  near  the  regions  of  interest  and  are  coarsely  distributed  elsewhere.  However,  these  non-uniform 
cells  can  sometimes  lead  to  grids  of  locally  poor  quality  that  will  increase  the  computational  solution  error.4 
The  performance  of  t  he  numerical  methods  can  be  significantly  reduced  even  if  just  a  few  low-quality  cells 
present.5 

The  grid  generation  efforts  can  be  traced  back  to  1960s*  Since  then  many  grid  generation  methods  have 
been  proposed  such  as  conformal  mapping,  algebraic  construction,  partial  differential  equation  solutions 
(elliptic,  hyperbolic,  and  parabolic  equations),  Delaunay,  advancing- front,  and  many  others.1  These  methods 
have  been  used  to  generate  structured  and/or  unstructured  meshes  over  simple  to  complex  geometries.  Each 
method  leads  to  a  different  grid  quality  to  the  same  geometry.  Besides,  while  some  mesh  generation  methods 
could  be  fast  and  even  automated,  others  are  a  time  consuming  manual  process. 

The  conformal  mapping  is  the  simplest  method  of  structured  grid  generation,  in  which  the  computational 
domain  is  mapped  onto  a  rectangular  region*  Although,  the  method  is  simple  and  efficient  and  creates  very 
high  quality  cells,  but.  it  is  practical  only  for  two-dimensional  simple  problems.  To  create  a  smooth  struc¬ 
tured  grid,  elliptic  grid  generators  can  be  used,  which  involve  the  numerical  solution  of  inhomogeneous  elliptic 
partial  differential  equations.6  Grid  smoothness  can  help  reducing  the  truncation  error  and  improving  the 
accuracy.  A  disadvantage  of  the  elliptic  grid  generation  is  the  limited  control  over  the  interior  grid  points. 
Algebraic  methods  have  also  been  used  to  create  structured  meshes*  These  methods  use  algebraic  transfor- 
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illation  and  an  interpolation  scheme  to  distribute  grid  points  from  a  discrete  set  of  data.  In  compai  ison  to 
elliptic  grid  generators,  algebraic  methods  require  much  less  computational  effort  and  has  bettei  control  over 
grid  point  locations,  but  the  grid  may  not  be  as  smooth  as  an  elliptic  grid,  io  improve  the  interior  mesh 
over  complex  geometries,  the  multi-block  strategy  has  been  proposed,  in  which  the  computational  domain 
is  divided  into  several  smaller  sub-domains  (block)  and  then  separate  meshes  are  generated  for  each  block.8 
Each  block  can  initially  be  meshed  by  an  algebraic  method  and  then  an  elliptic  solver  be  used  to  smooth 
the  meshes. 

Although,  structured  meshes  offer  higher  accuracy,  simplicity,  and  easy  data  access  compared  to  an 
unstructured  mesh,  but  unstructured  meshes  (using  tetrahedra  cells)  are  more  popular  for  complicated 
geometries.9  Delaunay  and  advancing- front  methods  are  the  most-  famous  tetrahedral  mesh  generators. 
The  advancing  front  creates  the  cells  one  by  one,  starting  from  the  domain  boundary  and  by  marching 
a  front  toward  the  interior.10  The  front  refers  to  the  cells  that  meet  unmeshed  domains;  the  front  will 
eventually  vanish  when  the  mesh  is  completed.  Advancing-front  approacli  will  have  the  best  quality  cells  at 
the  boundaries  and  the  worst  cells  where  the  front  collides  itself-  Therefore,  this  mesh  generator  is  veiy 
helpful  for  modeling  inviscid  or  laminar  flows  over  solid  walls.  The  advancing-front  method  can  even  be 
used  to  create  quadrilateral  cells,  however,  the  interior  cells  could  have  low  quality  or  the  mesh  cannot  be 
completed  for  a  complex  geometry.  Typical  Delaunay  mesh  generators  start  from  a  boundary  discretization. 
New  points  and  triangular  cells  are  then  added  to  satisfy  a  particular  connectivity.  I  lie  method  maximizes 
the  minimum  angle  of  all  triangles  to  avoid  sharp  and  distorted  cells  wherever  possible.11  Delaunay  mesh 
generators,  unlike  advancing-front  method,  has  the  worst  cells  at  the  domain  boundary  and  the  best  cells 
interior.10  Therefore,  the  Delaunay  and  the  advancing  front  approaches  are  often  combined  to  use  the 
advantages  of  both  methods. 

Unstructured  meshes  are  not  much  effective  when  localized  regions  of  high  gradients  appear  in  the  flow, 
such  as  boundary  lavers.12  Hybrid  grid  generators  and  anisotropic  tetrahedral  extrusion  are  commonly  used 
to  treat  boundary  layers  in  unstructured  meshes.  The  first  step  of  a  hybrid  approach  is  the  triangulation  of 
the  computational  domain  and  then  placing  prismatic  cells  near  boundary  surfaces.  The  prismatic  cells  are 
typically  created  using  an  advancing  layer  scheme.  In  such  a  mesh,  the  outside  boundary  nearly  has  isotropic 
tetrahedral  cells.13  A  disadvantage  is  that  prismatic  cells  may  have  poor  quality  near  sharp  edges.  An 
anisotropic  tetrahedral  extrusion  creates  triangular  cells  in  a  layer-by-layer  fashion  with  extrusion  direction 
normal  to  the  surface.12  The  method  begins  with  an  isotropic  tetrahedral  mesh  (generated  by  Delaunay 
and/or  advancing-front  methods),  and  tetrahedra  on  surfaces  are  then  subdivided  into  anisotropic  cells  for 
user-specified  number  of  layers.14  However,  this  approach  might  create  extremely  stretched  tetrahedra  cells 
near  the  Avails  which  impact  the  accuracy  of  computations. 

At  United  States  Air  Force  Academy  (USAFA),  the  hybrid  mesh  approach  has  been  successfully  used  for  a 
number  of  years  for  CFD  simulation  of  many  fighter  aircraft.15’16’17’18  These  grid  have  structured  cells  near 
solid  surfaces  and  unstructured  cells  elsewhere.  In  more  details,  an  inviscid  tetrahedral  mesh  is  generated 
using  the  ICEM-CFD  code.  This  mesh  will  then  be  used  as  a  background  mesh  by  the  mesh  generator  of 
TRITET19’20  which  builds  prism  layers  using  a  frontal  technique.  TRITET  rebuilds  the  inviscid  mesh  while 
respecting  the  size  of  the  original  inviscid  mesh  from  ICEM-CFD.  More  recently,  the  aut  hors  of  this  work  used 
this  hybrid  grid  approach  for  prediction  of  aerodynamic  characteristics  of  a  ram-air  parachute  21  ’ 22  The  grid 
sensitivity  study  showed  that  solutions  are  very  sensitive  to  the  grid  quality  and  size  for  the  air  foils/ wings 
with  an  open  inlet.  For  example  the  airfoil  solutions  became  grid  independent  for  the  hybrid  grids  that  have 
more  than  1.7  million  cells.  That  brings  to  mind  a  question  whether  another  grid  generator  could  produce 
similar  predictions  to  the  hybrid  fine  grids,  but  at  a  much  lower  computational  cost. 

The  objectives  of  the  work  are  twofold:  one  concerns  the  grid  resolution  and  to  find  coarse  grids  that  can 
have  comparable  results  to  the  fine  meshes.  Second,  is  to  understand  the  relationship  between  mesh  quality 
and  the  solution  accuracy  and  convergence.  Results  are  presented  for  two-dimensional  airfoils  with  open 
and  closed  inlet.  The  airfoil  sections  is  from  an  actual  parachute  canopy.  The  results  of  two-dimensional 
airfoils,  considered  in  this  work,  can  be  easily  generalized  to  three-dimensional  wings.  For  example,  a  high 
quality  3D  mesh  can  be  created  from  spanwise  extrusion  of  a  high  quality  2D  mesh.  This  work  is  organized 
as  follows.  First  the  Cobalt  flow  solver  is  described.  Several  grid  quality  metrics  are  also  reviewed.  Next, 
the  test  case  is  presented.  Finally,  the  simulation  results  will  be  discussed. 
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II.  CFD  Solver 


Cobalt  solves  the  unsteady,  three-dimensional,  compressible  Navier-Stokes  equations  in  an  inertial  ref¬ 
erence  frame.  Arbitrary  cell  types  in  two  or  three  dimensions  may  be  used;  a  single  grid  therefore  can  be 
composed  of  different  cell  types.23  In  Cobalt,  the  Navier-Stokes  equations  are  discretized  on  arbitrary  grid 
topologies  using  a  cell- centered  finite  volume  method.  Second-order  accuracy  in  space  is  achieved  using  the 
exact  Riemann  solver  of  Gottlieb  and  Groth,24  and  least  squares  gradient  calculations  using  QR  factoriza¬ 
tion.  To  accelerate  the  solution  of  the  discretized  system,  a  point-implicit  method  using  analytic  first-order 
in  viscid  and  viscous  Jacobians  is  used.  A  Newtonian  sub-iteration  method  is  used  to  improve  the  time  ac¬ 
curacy  of  the  point -implicit  method.  Toma.ro  et  al.2S  converted  the  code  from  explicit  to  implicit,  enabling 
Courant-  Fhedr  ichs-Lewy  (CFL)  numbers  as  high  as  10e.  In  Cobalt,  the  computational  grid  can  be  divided 
into  group  oi  cells,  or  zones,  lor  parallel  processing,  where  high  performance  and  scalability  can  be  achieved 
%*f!  011  fcf!‘  thousands  of  processors.26  Some  available  turbulence  models  in  Cobalt  are  the  Spalart-Allmaras 

^p?dd;“ '  Spahrt-Ailmaras  with  Rotation  Correction  (SARC),28  and  Delayed  Detached-eddy  simulation 
(DDES)  with  SARC.29 

Cobalt  checks  the  grid  quality  and  reports  a  score;  this  score  is  directly  related  to  a  particular  part  of 
t  le  second -older  accurate  spatial  operator  inside  Cobalt,30  The  reported  score  is  averaged  from  all  the  cells 
and  ranges  from  zero  to  hundred,  such  that  the  lower  the  grid  score,  the  more  numerical  dissipation  is  added 
to  the  solution  (worse  stability).  Note  that  Cobalt’s  grid  score  involves  the  mesh  geometry  information 
only,  and  not  tile  flow  solution  obtained  on  that  mesh.  Tiiis  means  that  the  grid  quality  is  a  fixed  number 
regardless  of  the  angle  of  attack  and  Mach  number.  Cobalt’s  User  Guide30  details  that  the  high  aspect  ratio 
of  cells  (typically  placed  in  the  boundary  layer)  can  cause  the  quality  to  suffer.  Also  regions  of  high  surface 
curvature  can  adversely  impact  grid  quality.  The  grid  quality  will  improve  if  a  sufficient  number  of  surface 
cells  are  used  to  accurately  capture  any  geometric  curvature. 


III.  Grid  Quality 


The  relationship  between  the  mesh  quality  and  solution  accuracy  will  be  investigated  in  this  work.  A 
piiori  grid  quality  metrics  could  provide  some  guidance  of  the  grid  before  running  it  in  CFD  33  The  mesh 
quality  mainly  deals  with  the  geometric  information  of  the  grid.  Before  examining  the  grid  quality,  one  should 
check  the  grid  for  negative  volume  (or  volumes  below  a  threshold)  and  folded  cells.  These  tvpe  of  cells  make 
a  cell-centered  flow  solver  impossible  to  iterate.32  Typically,  a  high  quality  mesh  is  obtained  by  creating  well- 
shaped  cells  (orthogonal  structured  cells  or  isotropic  tetrahedra  cells)  with  moderate  smoothness.  Therefore 
most  grid  quality  metrics  indicate  how  much  a  grid  cell  deviates  from  its  ideal  shape. 

A  great,  but,  more  advanced,  way  to  improve  the  grid  quality  metrics  is  to  add  information  about  the 
numerical  solution  obtained  on  the  mesh.4  One  simple  example  is  Y+  values  at  the  walls.  According  to  the 
gnc  c  ing  guidelines  from  the  2nd  AIAA  CFD  high  lift  prediction  workshop,33  approximate  initial  spacing 
normal  to  all  viscous  walls  should  have  Y+  values  of  approximately  1.0,  2/3,  4/9,  and  8/27  for  a  coarse! 
medium,  fine,  and  extra  fine  grid,  respectively.  In  more  advanced  grid  generation  methods,  the  solution 
errors  are  used  to  refine  the  mesh  in  regions  where  the  error  is  large. 

Mesh  resolution  and  quality  should  be  checked  prior  to  running  the  grid  in  CFD.  The  mesh  resolution 
should  be  high  enough  to  capture  the  flow  physics.  Mavriplis  et  a],34  for  example,  have  recommended  a 
chordwise  grid  spacing  of  0.1%  of  local  chord  (for  a  medium  grid)  at  the  wing  leading  and  trailing  edges  In 
the  span  wise  direction  at  the  wing  root  and  tip,  a  grid  spacing  of  0.1%  of  semispan  was  recommended  The 
number  of  recommended  cells  at  the  wing  trailing  edge  base  are  8,  12,  16,  24  for  a  coarse,  medium,  fine,  and 
extra  fine  grid,  respectively.  The  grid  resolution  can  be  evaluated  by  sensitivity  studies.  The  results  become 
mesh  independent  if  they  show  less  than  3%  difference  from  a  30%  finer  mesh.35 


The  grid  quality  is  often  measured  for  each  cell  from  given  information  about  the  cell’s  aspect  ratio  and 
skewness.  It  is  also  desirable  to  have  near- wall  faces  parallel  to  the  wall.  Besides,  the  rate  at  which  the  grid 
spacing  changes  from  one  cell  to  another  (grid  smoothness)  is  important.  Optimal  grids  have  equilateral 
cells  (equilateral  triangles  and  squares)  with  a  smooth  change  of  dimensions  through  the  domain.35  In  more 
details,  cell  skewness  should  be  kept  minimum.  High  aspect  ratios  would  slow  down  the  convergence.  Finally 
a  spacing  growth  rate  below  20%  does  not  affect  the  solutions,35  but,  higher  values  affect  the  solution  accuracy 
and  convergence.  The  reader  should  note  that  these  criteria  might  change  from  one  solver  to  another.  This 
is  discussed  iater  in  this  section. 
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There  are  some  available  grid  quality  metrics.  Alter,"*®  for  example,  has  described  a  giid  quality  metric 
for  structured  three-dimensional  grids.  This  metrics  is  defined  by, 


GQ  = 


(1) 


where  GQ  denote  the  grid  quality  ranging  from  0  to  one;  0milI.  0,  and  tmax  show  deviations  from  orthog¬ 
onality.  straightness,  and  stretching,  respectively.  Best  grid  quality  comes  with  GQ  =  i.  The  minimum  of 
the  average  skewness,  £>min  ,  is  determined  in  three  directions  by 


0  U— roll  slant-  t-OS 
7T 

2 

0  | constant  "T^OS 


|  C= constant^- 
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(3) 

(4) 


where  its  value  is  zero  and  minimum  for  a  collapsed  cell  and  is  one  and  maximum  for  an  orthogonal  edge 
segment.  The  straightness,  0,  measures  the  extra  distance  traversed  between  the  wall  and  the  outer  bound¬ 
ary  and  is  computed  normal  to  the  boundary  as 


0  = 


^  \  (Crtlin  Kmiix  ^ 


2^C=i 


(5) 


where  the  grid  line  straightness  ranges  from  zero  to  one,  where  zero  is  for  a  grid  line  that  begins  and  ends 
at  t  he  same  point.  Finally,  the  maximum  average  planar  stretching,  «illax,  measures  the  spacing  of  one  point 
to  the  next  in  all  three  directions  and  is  given  by 


max  (|rf  |+,ff  |_) 
min  ||jo|+,^|-) 
max(|f^|+.r^l~) 
min 

inax(|rc|  +  ,r<C) 


(6) 

(7) 

(8) 


where  each  term  ranges  from  one  as  minimum  to  infinity.  The  minimum  refers  to  a  uniform  spacing.  An 
abrupt  change  in  the  spacing  increases  the  stretching  value.  Alter’s  examples  of  orthogonality,  stretching, 
and  straightness  are  shown  in  Figure  1.  For  quadrilateral  (two-dimensional)  meshes,  Dannenh  offer31  defined 
below  quality  metrics: 
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.  min  (  di .  d2 ) 

skew  =  - — — — — 

max  (dj ,  d2) 

(9) 

,  -  ( s  1  ^3  s2  $ 4  \ 

taper  =  nun  — 

\  s2  *'4  J 

(10) 

aspect  ratio  -  +*3,*3 +  *) 

max  (&i  +  S3,  ^2  +  54) 

(11) 

stretch  -  JuOniAnAz^JU) 
v  max(Aj,  A2,  A3,.44) 

(12) 

wheie  si,  82,  S3,  and  s4  are  the  lengths  of  the  adjacent  sides  around  the  cell;  d.\  and  d2  are  the  lengths  of 
two  diagonals,  and  .4 1 ,  A2,  A3  and  .44  are  the  areas  of  four  cells  that,  surround  a  node.  All  these  metrics  are 
computed  at  eacli  cell  and  then  are  added  to  find  the  overall  grid  quality, 

A  grid  quality  metrics  based  on  the  deviation  from  an  ideal  cell  shape  will  probably  invalidate  most  of 
anisotropic  cells  near  the  walls  for  viscous  flows.  On  the  other  hand,  these  anisotropic  cells  are  allowed  and 
acceptable  in  many  flow  solvers.  In  addition,  for  all  cell-centered  flow  solvers,  the  quadrilateral  faces  should 
be  planar,  however,  non-planar  faces  are  not  an  issue  in  most  node-based  solvers.  McDaniel32  in  his  PhD 
dissertation,  therefore,  proposed  that  an  independent  mesh  quality  should  be  provided  from  the  viewpoint  of 
each  flow  solver  (cell- centered,  node-based  and  others).  He  then  introduced  grid  quality  metrics  for  the  kCFD 
code  '  which  solves  the  unsteady,  three-dimensional,  compressible  PANS  equations  on  hybrid  unstructured 
grids.  Note  that  both  Cobalt  and  kCFD  code  originated  from  the  Air  Vehicles  Unstructured  Solver  (AVUS. 
formally  known  as  CobaltbO24),  which  is  a  parallel,  implicit,  unstructured  flow  solver  developed  by  the  Air 
Force  Research  Laboratory.  Both  solvers  have  since  been  highly  modified,  McDaniel  then  defined  some  grid 
quality  metrics  which  are  applicable  to  the  most  cell-centered  numerical  algorithms  and  all  cell  types  and 
take  into  account  the  (act  that  large-area  faces  have  a  larger  flux  contribution  to  the  solution  value  in  the 
cell.  All  metrics,  also,  lie  in  the  range  from  zero  to  one,  while  one  represents  the  maximum  quality.  Some  of 
these  quality  metrics  are  briefly  reviewed  in  this  paper. 

Assume  that  there  are  no  folded  and  negative  volume  cells,  A  planarness  metric  is  defined  for  the  ith 
face  of  the  ith  cell  as 


l  planarj 


COS  1  (ri.j.1  j  1 
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face  nodes  <  3 
face  nodes  —  4 


(13) 


where  and  are  the  unit-normal  vectors  to  each  lace.  The  cell  planarness  metric  is  then  calculated 
as: 


No.  faces 

Q  planar  y,  .  ^  ^  V  planar, j  (Id) 

^3  3  j=l 


A/  denotes  jt!l  face  area  of  the  ith  cell.  Figure  2  (a)  gives  the  planarness  metric  values  for  some  example 
cells.  The  skew- smoothness  metric  for  the  jth  face  in  the  ith  cell  is  also  defined  as 


skew-smooth  tj 


min  (  |  n  |  ^rifiht  J  I ) 

. max  ( I  fioft . J  1 1  |  f'ri  ght  .j  I) 


2  (1  ^right,j  ] ) 


(15) 
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Pic[t  j  and  Prfght ,,  are  the  vectors  from  the  adjoining  left  and  right  cell  centroids,  respectively,  to  the  centroid 
of  the  jth  face  of  the  cell.  i>loft  j  and  t'nght.j  are  the  corresponding  unit  vectors.  The  first  and  second  terms 
of  equation  15  correspond  to  the  degree  of  smoothness  and  skewness  across  a  particular  face,  respectively. 
The  cell  skew-smoothness  metric  can  be  calculated  using  eq.  14.  If  the  cell  is  folded  then  y  skew- smooth  3 j  ^ ' 
The  skew-smoothness  metric  values  for  some  example  cells  are  given  in  Figure  2  (b), 

A  flow  alignment  (wall  orthogonality  or  straightness)  metric  was  also  defined  to  measure  the  degree  to 
which  cell  faces  are  parallel  to  t  he  solid  walls.  The  flow  alignment  metric  for  the  jth  face  in  the  i:l1  cell  is 


cos  1  (abs  ^hwaii-fi  face  j  )  ) 

Malign  j  =  _ 

j  aligns  =  2  [max  (ylligej  .  1  -  CtignJ )  —  0.5] 


(1C) 

(17) 


JIwaii 's  the  unit  normal  of  the  closest  solid  wall  boundary  face  and  Hfacej  is  is  the  unit  normal  of  the  j  face. 
Some  flow  alignment  metric  examples  are  shown  in  Figure  2  (c).  Other  quality  metrics  include  isotropy  and 
spacing  metrics  which  measure  the  deviation  from  an  “equilateral”  cell  shape  and  uniformity  in  the  mesh, 
respectively.  The  global  grid  quality  is  then  estimated  by 


No  of  metrics 

GQ  =  Pk 

k=l 


"No  of  cells 

y:  uj  (5jb,dwaiuj  Qk 

.  1=1 


(18) 


where  Pk  is  the  global  weight  for  the  kth  quality  metric;  </fe’  is  the  kth  quality  metric  value  for  the  ith 
cell;  w  ^st,dwaii.i)  is  a  weighting  function  for  each  quality  metric  which  depends  on  through  parameter  sk 

and  (iwaii  .■ ,  the  minimum  distance  of  the  i1*1  cell  centroid  from  the  solid  wall  boundary,  normalized  by  the 
characteristic  length  for  the  mesh. 


TV.  Test  Cases 

The  U.S.  DoD  program  to  develop  precision  guided  airdrop  systems  is  known  as  the  Joint  Precision 
Airdrop  System  (JPADS)  which  is  coordinated  by  the  U.S.  Army  Natick  Soldier  Center  (NSC).  The  JAPDS 
uses  round  and  ram-air  parachutes  (parafoils)  for  deceleration  and  control  of  the  payload.38  Different  type 
of  airfoil  sections  have  been  used  for  parafoils;  initial  ram-air  canopies  used  the  Clark-Y  section  which  has 
good  lift- to- drag  (L/D)  characteristics  for  medium  Reynolds  number,  Rt,  flows. 3 J  Over  the  years  this  airfoil 
shape  has  been  modified  for  parachute  applications  to  improve  lift-to-drag  ratio;  recent  canopy  designs  are 
based  on  the  airfoil  sections  used  in  glider  design  (for  example  NASA  LSI-0417  airfoil).40  In  this  work, 
low-speed  airfoil  section  of  an  actual  parafoil  is  investigated.  The  airfoil,  considered  here,  is  a  non- symmetric 
airfoil  that  has  a  flat  bottom  surface.  The  flow  around  this  airfoil  lias  been  studied  by  Ghoreyshi  et  al.2L22 
for  opening  and  closing  the  inlet;  the  open  and  dosed- inlet  geometries  are  shown  in  Figure  3. 

Ghoreyshi  et  al.21  performed  the  mesh-independent  study  of  both  airfoils.  These  results  are  shown  in 
Figures  4  and  5.  Closed  airfoil  grids  are  fully  structured  or  hybrid;  structured  grids  were  generated  using  the 
multi-block  techniques  inside  tire  commercial  code  of  Pointwise  V17.01.R3.  Hybrid  meshes  were  generated 
using  ICEM-CFD  code  and  the  mesh  generator  of  TRITET.  All  the  meshes  labeled  in  Figure  5  are  hybrid 
type. 

Figure  4  shows  that  solutions  of  structured  medium  (around  663,000  cells)  and  fine  meshes  (around  3 
million  cells)  match  everywhere.  The  solutions  ol  the  hybrid  medium  mesh  (around  744,000  cells)  are  in 
dose  agreement  with  the  medium  and  fine  structured  meshes  as  well.  The  coarse  structured  mesh  (around 
172.000  cells),  however,  only  matches  the  solutions  up  to  an  angle  of  6°;  at  higher  angles,  it  predicts  smaller 
lift  and  larger  drag  coefficients.  The  closed  airfoil  is  further  investigated  in  this  work.  This  geometry  is 
used  to  determine  if  predictions  of  a  low- resolution,  but  high-quality  mesh  can  still  match  with  the  medium 
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and  fine  grid  data.  The  second  goa!  is  to  relate  the  solution  accuracy  and  convergence  with  grid  parameters 
and  quality  metrics.  All  meshes  of  this  work  were  generated  using  the  Pointwise  code  and  predictions  were 
compared  with  solutions  of  the  fine  structured  grid  (containing  3M  cells).  Cobalt  reports  an  averaged  grid 
quality  of  99,74  in  a  range  of  zero  to  100  for  this  fine  grid, 

F lgure  5  shows  that  CFD  solutions  of  the  open-inlet  airfoil  largely  depend  on  the  grid  resolution.  Coarse 
and  medium  meshes  (containing  around  452,000  and  528.000  cells)  were  generated  from  a  low-density  inviscid 
grid  and  CFD  data  using  these  grids  do  not  match  with  fine  meshes  data  (containing  around  943,000.  1.7 
million,  and  2.45  million  cells).  Figure  5  shows  that  solutions  do  not  change  with  grid  density  for  the  grids 
that  have  more  than  1.7  million  cells  (see  linel  and  fine2  plots  in  the  figure).  The  open-inlet  airfoil  is  also 
used  in  this  work.  The  coarse  meshes,  with  different  grid  generation  methods,  will  be  tested  in  Cobalt,  and 
the  results  will  be  compared  to  finel  mesh  data  in  Ref.  21.  This  grid  has  a  Cobalt  quality  of  99.83.  The 
meshing  study  will  be  performed  on  both  outside  and  inside  domains.  The  effects  of  meshing  on  the  solution 
accuracy  and  convergence  will  be  investigated. 


V.  Results  and  Discussion 

For  the  sake  of  convenience,  the  closed  and  open  inlet  grids  are  labeled  -;CG"  and  i;OG”,  respectively 
All  CFD  simulations  were  performed  using  the  Spalart-Allmaras  (SA)  turbulence  model  and  were  run  on  the 
Cray  XE6  machine  at  the  Engineering  Research  Development  Center  (ERDC).  The  free-stream  velocity  in  all 
simulations  was  fixed  at  Maclr  0.1  and  the  Reynolds  number  is  1.4  xlO0  at  standard  sea  level  conditions.  The 
simulation  were  performed  for  an  angle-of-attack  sweep  from  zero  to  ten  degrees  with  one  degree  increment. 
In  all  simulations,  second-order  accuracy  in  time,  3  Newton  sub- iterations,  and  20,000  iterations  with  a  time 

step  of  le-4  seconds  were  used.  The  last  10,000  iteration  values  were  averaged  to  obtain  overall  lift,  and  drag 
coefficients. 

All  grids  were  run  in  Cobalt  and  errors  in  the  force  coefficients  from  the  predictions  of  fine  meshes  are 
estimated.  In  more  details,  the  error  norm  of  lift  coefficient  (CL),  drag  coefficient  (CD),  and  lift  to  drag 
ratio  (L/D)  is  defined  as: 


V^Ef=i(<cwGlid  -  yJi!lcGrid)2 

|2/FincGrid(moa.)   yFincGrid(mjn) 


x  100 


(19) 


where  y  =  [C'lX'd-L/D]:  N<  is  number  of  angles  of  attack  which  is  11  in  this  work.  The  error  norm  is  also 
found  for  the  maximum  lift  coefficient  (CLtrmK}  as 


err  — 


n,  NcwGrid  _  io  FincGrid 

FincGrid 

^  Lmax 


x  100 


(20) 


First,  results  are  presented  and  discussed  for  the  closed-inlet  grids.  A  medium  size  grid  (named  CGI) 
was  generated  with  details  given  in  Table  1.  The  grid  was  generated  layer  by  layer  using  an  elliptic  extrusion 
method  starting  from  the  walls.  The  method  stops  when  the  new  grid  layer  reaches  a  total  height  of  200m. 
The  outer  edges  of  the  cells,  at  the  final  layer,  define  the  freest  ream  boundary.  The  coarse  grid  used  in 
Rel.  21  was  also  modified  to  Stave  exactly  the  same  number  of  grid  points  and  spacings  on  the  wails  as  the 
CGI  mesh  has.  This  new  grid  (named  CG2),  however,  was  generated  using  a  multi-block  method  with  an 
algebraic  grid  generator  applied  to  each  block.  Grid  details  are  also  given  in  Table  1.  Table  1  confirms  that 

CGI  and  CG2  grids  have  nearly  equal  number  of  cells.  The  grid  overview  near  trailing  edge  is  shown  in 
Figure  6. 

Table  2  compares  Y+,  Cobalt  grid  quality,  and  the  error  norms  of  CGI  and  CG2  grids.  While  both  grids 
have  approximately  the  same  resolution,  Y+,  and  grid  quality,  but  the  errors  from  CGl  are  one  order  of 
magnitude  less  than  CG2  grid.  In  more  details,  Figure  7  compares  the  lift  and  drag  coefficients  of  CGI  and 
CG2  grids  with  the  fine  grid  data  of  Figure  4.  The  comparisons  show  that  CGI  grid  perfectly  matches  with 
the  fine  grid  predictions,  however,  CG2  predictions  do  not  match  everywhere.  A  question  that  arises  is  why 
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these  grids,  both  structured,  with  the  same  7+,  grid  quality,  and  resolution  lead  to  different  answers  to  the 
exact  same  problem. 

Prior  to  answer  this  question,  let's  first  examine  the  convergence  solution  of  CG 1  and  CG2  grids,  t  igure  5 
shows  the  convergence  histories  oi'  the  density  residual  {Dp/Dt)  for  CGI  and  CG2  grids  as  well  as  the  fine 
mesh  containing  3M  cells.  Prom  this  figure,  it  can  be  seen  that  both  CGI  and  CG2  grids  reached  converged 
values  with  approximately  3.0  order  of  magnitude  reduction  in  the  density  residual.  Figure  8  also  indicates 
that  both  grids  exhibit  very  similar  convergence  behavior.  Since  CGI  and  CG2  have  the  same  quality  and 
resolution,  a  preliminary  conclusion  is  that  solution  convergence  in  Cobalt  is  related  to  the  grid  quality  and 
the  resolution;  this  will  be  further  examined  later  for  more  grids.  The  fine  mesh  quality  in  Cobalt  is  99.74 
slightly  higher  than  CGI  and  CG2;  also  it  is  a  very  high  spatial  resolution  grid.  As  expected,  Figure  8  shows 
that  the  fine  mesh  has  better  convergence  than  other  grids. 

Now,  in  attempting  to  answer  the  aforementioned  question,  the  reader  is  referred  again  to  the  mesh 
overviews  shown  in  Figure  6.  Detailed  visualization  of  the  cells  around  the  trailing  edge  reveals  that  CGI 
has  better  wall  orthogonality  (straightness),  smoothness,  and  skewness  compared  to  CG2.  Many  cells  m 
CG2  are  highly  skewed;  besides  most  mesh  lines  are  not  orthogonal  to  the  walls.  However,  Cobalt  s  overall 
grid  qualitv  is ‘very  similar  for  both  grids.  The  grid  quality  plots  of  CGI  and  CG2  (around  the  trailing  edge) 
can  also  be  seen  in  Figures  9  (a)  and  (b).  Note  that  white  color  cells  in  the  figures  indicate  high  grid  quality 
(above  99).  Figure  9  (b)  shows  that  most  CG2  cells  in  Cobalt  have  high  quality  while  they  were  expected 
to  have  low  quality  because  of  the  poor  flow  alignment  and  skewness. 

Grid  smoothness  and  skewness  values  from  Pointwise  are  also  shown  in  Figure  9.  This  figure  shows 
that  CGI  cells  are  well-shaped  with  good  smoothness.  Comparing  Pointwise  plots  with  Cobalt  grid  quality 
pictures  reveals  that  Cobalt  grid  quality  is  related  to  the  grid  smoothness,  but  not  the  skewness  or  the  wall 
orthogonality.  Figure  10  shows  the  same  correspondence  for  the  cells  near  the  leading  edge.  CG2  cells  at 
the  leading  edge  have  slightly  better  smoothness  than  CGI,  but  it  has  highly  skewed  cells.  Figure  10  shows 
that  Cobalt’s  quality  (at  the  leading  edge)  is  better  for  CG2  compared  to  C’Gl  cells.  This  again  confirms 
that  Cobalt  grid  quality  does  not  change  with  the  skewness  or  wall  orthogonality.  This  means  that  a  high 
grid  quality  in  Cobalt  corresponds  to  good  smoothness  which  would  improve  the  convergence  but  it  does  not 
necessarily  provide  a  better  accuracy. 

Four  new.  grids  (CG10  -  CG13)  were  generated  around  the  closed-inlet  airfoil  by  the  normal  extrusion  ol 
the  wall  using  an  elliptic  mesh  generator.  The  edge  lengths  oi  these  grids  are  about  half-size  of  the  edges  in 
CGI  and  CG2  grids.  All  grids  are  structured  and  contain  around  120,000  cells.  The  near- wall  grid  spacing 
(Asl)  is  4e-5m  for  these  grids  which  makes  overall  V+  near  one.  The  difference  between  these  four  grid  is 
only  due  to  the  farfield  length  (/);  it  varies  from  25  to  100  chord  lengths.  Grid  details  and  errors  (with  the 
fine  mesh)  are  given  in  Tables  1  and  2. 

Figure  11  (a)  shows  that  CG10-CG13  grids  have  very  similar  convergence  behavior:  again,  because  these 
grids  have  the  same  resolution  and  quality.  Figure  11  (a)  shows  that  the  farfield  length  had  no  considerable 
effect  on  the  convergence  rate  in  Cobalt.  Figure  11  (b)  also  shows  the  error  trends  with  farfield  length  for 
grids  CG10  -  CG13.  Increasing  the  farfield  length  above  50  chord  length  does  not  change  much  the  Cl, 
Co  and  CLlnax  errors;  Figure  11  (b)  shows  that  L/D  predictions  can  be  improved  by  making  the  farfield 
boundary  bigger  but  its  impact  becomes  smaller  above  50c.  For  all  subsequent  grids,  a  farfield  length  of  50 
chord  will  therefore  be  used. 

CG2R-CG21  grids  were  generated  by  similar  mesh  generation  methods,  however,  the  near- wall  grid  spacing 
(Asi)  varies  in  these  grids.  The  convergence  data  of  these  grids  (CG20  -  CG21)  are  shown  in  Figure  12  (a) 
which  shows  the  convergence  improves  as  decreases  from  one  to  around  0.2  (CG11  in  this  figure  has 

=  1)+  The  solution  accuracy  is  also  improved  with  decreasing  Y*  as  shown  in  12  (b).  Figure  13  also 
compares  the  lift  and  drag  coefficient  predictions  of  CG21  grid  =0.2663  and  containing  134,000  cells) 
with  predictions  from  CGI  (T+=0-2663  and  containing  34LOOO  cells).  The  results  show  that  CFD  data  of 
the  coarse  grid  match  very  well  with  the  medium  grid  data.  Note  that  CGI  data  also  match  with  CFD  data 
of  the  fine  grid.  These  results  show  that  a  grid  size  of  around  134,000  cells  and  V  +  —  0,26  with  well-shaped 
cells  will  match  the  fine  data  of  Ref.  21. 

Three  unstructured  meshes  (CG30-32)  were  also  considered;  these  meshes  w^ere  generated  by  the  Delaunay 
I  et railed ra liz at io n  and  have  anisotropic  cells  near  the  wall  with  a  growth  rate  of  1.1.  These  grids  are  very 
similar,  except  that  they  have  different  initial  spacing  near  the  wall.  Gird  details  and  errors  are  given  in 
Tables  1  and  2.  These  grids  have  around  115,000  cells  and  have  grid  quality  above  99,  Figure  14  compares 
CG21  (structured)  with  the  CG32  (unstructured).  The  anisotropic  cells  near  the  wall  can  be  seen  in  this 
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figure.  Note  that  anisotropic  tetrahedral  layers  are  not  constant  everywhere;  the  local  extrusion  at  each 
point  would  stop  if  the  cell  size  become  dose  to  the  size  of  outside  isotropic  cells. 

I  iguie  15  compares  cobalt  quality  plot  of  CG32  with  the  Pointwise  skewness  quality  picture  of  CG32. 
figure  15(b)  shows  that  the  cells  near  the  wail  are  extremely  skewed,  but  Cobalt  considers  these  cells  high- 
quality  ones  because  of  the  grid  good  smoothness.  The  convergence  history  of  CGI,  CG21.  and  CG32  are 
plotted  in  Figure  16.  Despite  having  a  poor  skewness  near  the  wall,  CG32  has  very  similar  convergence 
behauor  to  CG_1  which  is  a  structured  grid.  This  again  confirms  that  Cobalt  convergence  does  not  change 
much  with  the  skewness  quality.  CGl  has  better  convergence  than  the  other  grids  because  it  is  a  high- 
resolution  compared  with  CG21  and  CG32. 

Although  CG32  and  CG21  have  similar  quality  and  convergence,  but  Table  2  shows  that  CG32  has  much 
arger  eiiors  m  the  li ft.  and  drag  coefficients  than  CG21.  In  more  details.  Figure  17  compares  the  lift  and 
drag  values  for  grids  of  CGI,  CG21,  and  CG32.  Figure  17  shows  that  CG32  match  well  with  CG21  and 

CG™o.  a!'  SmaU  ?ngles  of  attack;  however,  for  angles  above  6°,  CG32  lift  overestimates  CG21/CG1  data 
and  -G3_  drag  underestimates  them.  Finally,  for  the  structured  grids  of  CGll,  CG20,  and  CG21.  Table  2 
shows  that  the  accuracy  is  improved  by  decreasing  Y+  from  one  to  0.26;  this,  however,  does  not  apply  to 
the  unstructured  grids.  Decreasing  the  initial  spacing  makes  the  cells  near  the  wall  more  skewed  assuming 
the  wall  spacing  is  unchanged.  Poor  skewness  would  impact  the  solution  accuracy. 

Final  closed  girds  considered  are  hybrid  type  with  prismatic  layers  near  the  wall  and  isotropic  cells 
elsewhere.  Three  hybrid  grids  of  CG40,  CG41,  and  CG42  were  generated.  The  prismatic  layers  were 
generated  by  the  wall  normal  extrusion  using  an  elliptic  solver.  The  outer  domain  is  then  meshed  by  the 
Delaunay  tetrahedralization.  These  grids  are  much  coarser  than  the  structured  grids  of  CGI  0  to  CG21  Gird 
details  and  errors  are  gain  given  in  Tables  1  and  2.  The  main  difference  between  these  grids  is  the  number 
of  prism  layers.  Figure  18  shows  the  mesh  overview  of  CG2  with  75  prism  layers.  The  effects  of  the  number 
prism  layers  on  the  errors  can  be  seen  in  Figure  19  which  shows  the  accuracy  is  significantly  improved  by 
increasing  the  layers  from  25  to  50:  however,  the  accuracy  does  not  change  much  by  increasing  layers  from 
rr,°  ,  S  ^“fmpares  the  Iift  and  draS  coefficients  of  the  hybrid  CG42  mesh  with  predictions  of 
I.G1  and  CG21  grids.  The  results  show  that  predictions  from  a  hybrid  grid  with  75  prism  layer  still  match 
quite  well  with  expected  data. 

Next,  results  are  presented  for  the  open-inlet  geometry.  The  computational  domain  of  this  geometry 
consists  of  two  parts  corresponding  to  the  outside  and  inside  of  the  airfoil.  These  domains  are  meshed 
separately  in  this  work.  First  grids  considered  are  structured;  outside  mesh  was  generated  by  the  wall 
normal  extrusion  and  the  inside  mesh  is  an  algebraic  mesh  nearly  uniform  in  spacing  even  at 'the  inside 

M  °G5'  AU  these  «***  have  the  ^side  mesh  containing  around 

_r8,000  cel  Is  inside  OGl-  OG3  have  the  initial  spacing  of  4e-5m  at  the  outside  walls  but  they  have  different 
num  jer  of  layers  of  constant  spacing.  OG4  and  OG5  have  5  layers  of  constant  spacing  at  the  outside  walls 
but  the  initial  spacing  of  these  meshes  are  2e-5m  and  le-5m,  respectively.  OGl  mesh  at  the  leading  and 
trailing  edges  of  the  airfoil  is  shown  in  Figure  21.  OGl  has  25  layers  of  constant  spacing  at  the  outside  walls 
which  can  easily  be  seen  in  Figure  21 

nr-?G1  gnCfS  contam  around  400,000  cells.  Figure  22  (a)  shows  the  convergence  histories  of  OGl 

QG_,  and  OG3  grids  compared  with  the  hybrid  fine  mesh  (1.7M  cells)  convergence.  The  convergence 
comparisons  show  that  all  coarse  grids  converged  to  the  same  values;  again  due  to  having  similar  resolution 
am  quality.  In  comparison  to  the  closed-inlet  airfoil  convergence  plots.  Figure  22  (a)  shows  that  the  open- 
mlet  take  a  longer  time  to  converge.  Also,  Figure  22(a)  shows  that  the  fine  grid  convergence  is  not  significantly 
different  from  the  coarse  grids, 

,,,  effect®  of  *h€  mlmber  of  constant-spacing  layers,  at  the  wall,  (Ns )  on  the  errors  are  shown  in  Figure  22 
(  ).  he  results  show  that  Ns  has  small  effects  on  the  open-inlet  solutions.  In  more  details,  the  lift  and  drag 
coefficients  of  OGl ,  OG2,  and  OG3  grids  are  compared  with  the  fine  mesh  predictions  in  Figure  23.  The  lift 
and  drag  coefficients  are  split  into  the  inner  and  outer  surfaces.  Notice  that  inner  lift  does  not  change  with 
the  angle  of  attack:  Gboreyshi  et  al  1  showed  that  flow  is  stationary  inside  with  stagnation  pressures  almost 
everywhere.  Figure  23  (b)  however,  shows  that  the  outer  surface  has  negative  drags  before  stall  angle  and 
the  inner  surface  has  very  large  drag.  The  comparison  plots  of  Figure  23  show  that  coarse  grids  with  uniform 
spacing  inside  match  very  well  with  the  fine  mesh  force  coefficients  at  inside  and  outside  walls. 

OG2,  OG4,  and  OG5  have  around  278,000  cells  inside  and  very  similar  grids  outside;  in  these  grids 
V  ”  0  ^,ld  Sr»wth  rate  is  1.1.  However,  the  initial  grid  spacing  is  different;  Y+  values  of  these  grids  are 
given  in  Table  4  and  ranges  from  one  for  the  OG2  to  0.23  for  the  OG5.  The  effects  of  F+  on  the  errors  are 
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shown  in  Figure  24  and  compared  with  trends  of  closed-inlet  solutions.  Figure  24  shows  that  the  errors  of 
the  open-inlet  grids  are  slightly  reduced  by  decreasing  Y+ ,  however,  the  errors  oi  all  )  values  considered 
are  small  and  less  than  2%.  On  the  other  hand,  the  solution  accuracy  of  the  closed-inlet  grid  significantly 
changes  with  Y+. 

OG5,  OGIO,  OG11,  OG12  grids  were  used  to  study  the  effect  of  inside  mesh  resolution  on  the  overall 
solutions.  All  these  grids  nearly  have  very  similar  mesh  outside  but  inside  mesh  changes  from  278,000  to 
22.000  cells.  This  is  achieved  by  reducing  the  inside  wall  grid  dimensions.  OG12  has  only  22.000  cells  inside 
and  tire  gird  is  shown  against  OG5  with  278,000  cells  in  Figure  25.  Figure  26(a)  shows  that  OG12  will  less 
cells  has  density  residual  to  a  larger  value  than  OG5  and  the  fine  mesh;  however,  the  accuracy  is  still  as 
good  as  others  as  shown  in  Figure  26(b)  which  shows  OG12  errors  from  the  fine  mesh  are  less  than  1%.  This 
means  that  open-inlet  solution  is  riot  much  sensitive  to  the  number  of  cells  inside. 

Figure  27  shows  the  flow  solutions  of  the  closed  and  open  inlet  at  a  =  8.5°.  This  figure  shows  that 
the  flow  inside  airfoil  is  stationary  almost  everywhere;  therefore  the  flow  solution  should  not  change  much 
with  inside  mesh  resolution.  These  investigate  the  inside  mesh  type  on  the  solution  OG20  and  OG21  grids 
were  generated:  these  grids  again  have  structured  mesh  outside  same  as  previous  meshes,  but  the  inside  was 
meshed  by  the  Delaunay  tetrahedralization.  OG20  has  only  9.800  tetrahedra  cells  inside;  the  grid  is  shown 
in  Figure  28(a).  Table  4  and  Figure  29  show  that  QG20  with  coarse  tetrahedra  cells  inside  still  match  with 
the  fine  data.  Therefore  solutions  not  much  depend  on  the  inside  cell  resolution  as  well  as  type. 

OG22  mesh  has  unstructured  cells  inside  and  outside.  The  outside  is  meshed  by  the  Delaunay  tetrahe¬ 
dralization  and  has  anisotropic  layers  at  the  wall.  Inside  cells  are  uniform  isotropic  tetrahedra  cells.  The 
mesh  and  solutions  are  shown  in  Figure  28(h)  and  Figure  29,  respectively.  Figure  29  shows  that  outer  surface 
and  therefore  total  lift  and  drag  do  not  match  with  fine  data  everywhere.  This  is  probably  due  to  having 
highly  skewed  cells  near  the  outer  walls. 

Filial  grid.  OG40,  is  a  hybrid  mesh  outside  and  unstructured  inside.  Viscous  layers  at  outside  wall  were 
generated  bv  the  wall  normal  extrusion;  it  has  75  layers.  The  rest  of  outside  domain  is  meshed  by  the 
Delaunay  mesh  generator.  The  grid  overview  is  shown  in  figure  30.  The  lift  and  drag  coefficients  of  this 
grid  are  compared  with  the  fine  and  OGl  mesh  in  Figure  31  which  shows  hybrid  grid  match  well  with  other 
grids  data;  there  are  very  small  discrepancies  in  lift  coefficients  in  the  post-stall  region. 


VI.  Conclusions 

This  paper  provides  an  overview  on  the  grid  quality  and  resolution  effects  on  tlie  aerodynamic  modeling 
of  parachute  canopies.  The  CFD  simulations  were  performed  using  the  Cobalt  flow  solver  on  two  dimensional 
canopy  sections  with  open  and  closed  inlets.  The  simulation  results  show  that  Cobalt  grid  quality  (ranging 
from  zero  to  100)  depends  mainly  on  the  mesh  smoothness,  but  not  the  skewness  or  the  wall  orthogonality. 
Cobalt  convergence  is  improved  with  increasing  grid  resolution,  grid  quality  (smoothness)  and  decreasing 
initial  grid  spacing  at  the  wall  (Y 4  ).  The  results  showed  that  while  a  high  quality  grid  of  Cobalt  leads  to 
a  better  convergence,  but  is  does  not  always  lead  to  a  better  accuracy.  Ilie  solution  accuracy  will  change 
with  grid  resolution  and  smoothness,  as  well  as  skewness  and  the  wall  orthogonality. 

The  results  of  this  work  showed  that  the  solutions  of  coarse  structured  and  hybrid  grids  (around  100,000) 
match  well  with  prediction  of  fine  meshes  containing  a  few  million  cells,  i'hese  grids  have  appropriate  Y+ 
values,  good  smoothness,  skewness,  and  wall  orthogonality.  While  unstructured  meshes  with  anisotropic  cells 
near  the  wall  have  very  good  grid  quality  in  Cobalt  but  they  have  the  worst,  accuracy  between  considered 
grids  because  of  poor  skewness  at  the  walls.  The  results  also  showed  that  in  comparison  to  the  closed  inlets, 
the  open  geometry  solutions  are  less  sensitive  to  the  initial  grid  spacing  and  number  of  constant  spacing 
layer  at  outside  walls.  Also,  the  solutions  do  not  change  with  inside  mesh  resolution  and  type.  However,  a 
mesh  with  anisotropic  layers  at  outside  walls  have  again  the  worst  accuracy  between  considered  grids. 

The  results  of  two-dimensional  airfoils,  considered  in  this  work,  can  be  easily  generalized  to  three- 
dimensional  wings.  For  example,  a  high  quality  3D  mesh  can  be  created  from  spanwise  extrusion  of  a 
high  quality  2D  mesh.  Our  future  work  will  expand  the  results  to  include  meshes  around  three-dimensional 
wings  with  open  inlet,  deflected  trailing  edges,  and  bleed-air  spoilers.  The  experiments  of  these  wings  are 
underway  at  the  USAFA  subsonic  wind  tunnel. 
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I  able  1.  Closed  gnd  details.  A,i  and  N,  denote  the  initial  spacing  and  the  number  of  layers  of  constant 
spacing  normal  to  all  viscous  walls,  respectively.  Growth  rate  in  the  viscous  layer  is  shown  by  GR.  Finally,  c 
is  the  airfoil  chord  and  l  denotes  the  farheld  length  away  from  the  airfoil* 


Grid 

Type 

Method 

A„i 

Ns 

GR 

l/c 

#  Cells 

CGI 

STR 

Normal  extrusion t  elliptic 

lc-5  m 

7 

1.1 

100 

341,850 

CG2 

Sl'R 

Multi- block,  algebraic 

lc-5  m 

7 

1.1 

100 

324,570 

CG10 

STR 

Normal  extrusion,  elliptic 

4e-5  m 

- 

LI 

25 

1 15,320 

CGll 

STR 

Normal  extrusion,  elliptic 

4e-5  m 

LI 

50 

120,900 

CG12 

STR 

Normal  extrusion,  elliptic 

4e-5  m 

_ 

1.1 

75 

124,620 

CG13 

STR 

Normal  extrusion,  elliptic 

4e-5  m 

- 

171 

100 

128,340 

CG20 

STR. 

Normal  extrusion,  elliptic 

2e-S  in 

_ 

LI 

50 

127,410 

CG21 

STR 

Normal  extrusion,  elliptic 

le-5  m 

. 

LI 

50 

133,920 

CG30 

UN  STR 

Delaunay,  anisotropic  viscous  layers 

4e-5  m 

_ 

Id 

50 

101,894 

CG31 

UN  STR. 

Delaunay,  anisotropic  viscous  layers 

2c-5  m 

_ 

LI 

50 

114,970 

CG32 

UNSTR 

Delaunay,  anisotropic  viscous  layers 

le-5  m 

. 

LI 

50 

128,564 

CG40 

HYBRID 

Delaunay,  25  prismatic  viscous  layers 

lc-5  ni 

_ 

Id 

50 

55,324 

CG41 

HYBRID 

Delaunay,  50  prismatic  viscous  layers 

le-5  m 

. 

Id 

50 

77,754 

CG42 

HYBRID 

Delaunay,  75  prismatic  viscous  layers 

Lc-5  m 

- 

LI 

50 

100 .84 2 

Table  2-  Closed  grid  solution  details* 


err  % 

Grid 

r+ 

Quality 

CL 

CD 

LjD 

0^711  Q  3T 

CGI 

0.2663 

98,62 

0,1484 

0.1276 

0.2209 

0.1183 

CG2 

0,2339 

98.81 

1.1504 

2.0752 

5.6311 

-1.7914 

CG10 

1,074 

99.14 

1.8921 

2.9175 

8  1312 

-1,7116 

CGll 

1,074 

99,15 

L7397 

2.3424 

5,0274 

-L6729 

CG12 

1.074 

99, 16 

1.7155 

2.2623 

4.4330 

-L6793 

CG13 

1.074 

99,19 

1.7045 

2.2196 

4.0466 

-1,6824 

CG20 

0.5336 

99,12 

L0534 

1,1167 

2,7766 

-1,0541 

CG21 

0,2663 

98,99 

0.8772 

0.7775 

2.2576 

-0.8549 

CG3Q 

0.7260 

99,83 

0,3247 

1  *8036 

4.3489 

0*3532 

CG31 

0.3642 

99.71 

4*3107 

5.8003 

3.6808 

1.4514 

CG32 

0.1796 

99,55 

3.5703 

4,5855 

2,9438 

1.6056 

CG40 

0,2663 

98.55 

5.2646 

7,2829 

7.8631 

-6,2807 

CG41 

0.2663 

98.97 

1.0484 

0.5219 

1.6164 

-0.8507 

CG42 

0.2663 

98.94 

0.5906 

0,8316 

2dl28 

-0.5507 

14  of  29 


American  Institute  of  Aeronautics  and  Astronautics 


Table  3.  Open  grid  details,  ijc  =100  and  GR  —  1,1  for  all  grids. 


Grid 

Outside  mcsli 

Inside  mesh 

Afii 

N* 

#  Cavity  colls 

#  Total  cells 

0G1 

STR L,  elliptic 

STR,  uniform 

4g-5  m 

25 

278,520 

421,840 

OG2 

STR,  elliptic 

STR,  uniform 

4e-5  m 

5 

278,620 

403.240 

OG3 

STR,  elliptic 

STR,  uniform 

4  c -5  m 

0 

278,620 

399,520 

OG4 

STR.,  elliptic 

STRt  uniform 

2e-5  m 

5 

278,620 

407,500 

OG5 

STR ,  elliptic 

STR,  uniform 

le-5  m 

5 

278,620 

416/260 

OGIO 

STR,  elliptic 

STR,  uniform 

lc-5  m 

5 

120,384 

285,028 

OG11 

STR ,  elliptic 

STR,  uniform 

le-5  m 

5 

56,000 

194,188 

OG 12 

STR,  elliptic 

STR,  uniform 

le-5  m 

5 

22,000 

148,262 

OG20 

STR,  elliptic 

UN  STR,  Delaunay 

le-5  m 

5 

9,812 

137,058 

OG21 

STR,  elliptic 

UN  STR,  Delaunay 

lc-5  m 

5 

5.209 

136,520 

OG22 

isotropic  and  anisotropic 

UN  STR,  Delaunay 

le-5  rn 

- 

5,275 

118,202 

OG40 

HYBRID,  75  prism  layers 

UNSTR,  Delaunay 

le-5  m 

- 

5,275 

98,620 

Table  4.  Open  grid  solution  details. 


err  % 

Grid 

Y+ 

Quality 

CL 

Co 

L/D 

Cxm  O.T. 

OG  1 

0.9490 

99.66 

0.7146 

1,1728 

1,7029 

-0.4407 

OG2 

0.9495 

99.69 

0.9527 

1.4604 

1,9284 

-0.5427 

OG3 

0.9501 

99.70 

1,2356 

1.6935 

2,3153 

-0.7218 

OG4 

0,4721 

99.66 

0.4549 

0.5146 

0.5746 

-0.1788 

OG5 

0.2360 

99,59 

0.4232 

0.3767 

0.4536 

-0,1312 

OGIO 

0.2360 

99.26 

0,4269 

0.3506 

0.4322 

-0.0934 

OG11 

0,2359 

98.85 

0,4430 

0.4462 

0.5599 

-0.1770 

OG  12 

0.2357 

98.86 

0.5402 

0,4694 

0.S619 

-0,4866 

OG20 

0.2360 

98,92 

0,4994 

1.5627 

2.6723 

-0.1016 

QG21 

0.2360 

98.91 

0,4323 

1.8874 

3,8378 

-0.1263 

OG22 

0.1599 

99.35 

3.4452 

1,1538 

2.6855 

1.7791 

OG40 

0.2364 

98.90 

1,0078 

1.5233 

3,2509 

0.3833 

— _ 
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(a)  Poor  orthogonality  (b)  Good  orthogonality 


(e)  Poor  straightness 


(f)  Good  straightness 


Figure  1.  Orthogonality,  stretching,  and  straightness  examples  for  structured  grids. 


36 
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q  =  1 .0 


q  =  0.80 


q  =  0,50 


q  =  &.  1 1 


(a)  planarness  quality  metric  examples 


(b)  skew-smoothness  metric  examples 


e£ZCZ3 


q  =  0.52  q-0"3  q~UQ 

(c)  flow  alignment  metric  examples 


Figure  2.  Planarness,  skew -smoothness,  and  flow  alignment  metric  examples  from  Ref.  32  * 
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y(m) 


(a)  Closed  inlet  (b)  Open  inlet 

Figure  3*  Airfoil  geometries  for  the  parachute  cases 


Figure  4.  Mesh-i independent  study  of  the  closed- inlet  airfoil. 


(a)  Lift,  coefficient  (b)  Drag  coefficient 


Figure  5.  Mesh-independent  study  of  the  open-inlet  airfoil. 
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(a)  CGI,  structured,  elliptic  mesh  (b)  CG2.  multi-block  structured  mesh 

Figure  6.  CGI  and  CG2  grids  overview* 


Figure  7.  Lift  and  drag  coefficients  of  CGI  and  CG2  grids.  The  structured  fine  grid  is  from  Ref.  21 


Figure  8.  Density  residual  convergence  (Dp/Dt)  of  CGI  and  CG2  grids,  where  p  and  t  denote  density  and 
time,  respectively*  The  structured  fine  grid  is  from  Ref*  21. 
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(a)  CGI,  Cobalt  grid  quality 


(b)  CG2,  Cobalt  grid  quality 


(c)  CGI,  Pointwisc  smoothness 


(d)  CG2,  Pointwisc  smoothness 


(c)  CGI,  Pointwisc  skewness 


(f)  CG2,  Pointwisc  skewness 


Figure  9.  Cobalt  grid  quality  relationship  with  mesh 
the  white  colored  cells  have  the  highest  quality. 


geometry-  Mesh  is  around  trailing  edge. 


In  (a)  and  (b) 
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(a)  CGI,  Cobalt,  grid  quality 


(c)  CGI,  Pointwisc  smoothness 


(c)  CGI,  Pointwisc  skewness 


(b)  CG2,  Cobalt  grid  quality 


(d)  CG2h  Pointwisc  smoothness 


(f )  CG2,  Pointwisc  skewness 


Figure  10,  Cobalt  grid  quality  relationship  with  mesh  geometry-  Mesh  is  around  leading  edge.  In  (a)  and  (b) 
the  white  colored  cells  have  the  highest  quality. 
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(a)  density  residuals 


Figure  11-  Effects  of  far  field  length  (l)  on  the  solution  convergence  and  errors, 


(a)  density  residuals 


Figure  12,  Effects  of  initial  grid  spacing  (Asi)  on  the  solution  convergence  and  errors  of  structured  grids. 


Figure  13.  Lift,  and  drag  coefficients  of  CGI  (341,850  cells)  and  CG21  (133,920  cells)  grids. 
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(a)  CG21,  STR.,  elliptic  mesh 


(b)  CG32.  UN  STR  with  anisotropic  viscous  layers 


Figure  14,  CG21  and  CG32  grids  overview. 


(a)  Cobalt  grid  quality  (b)  Mesh  skewnewss 

Figure  15.  Cobalt  cell  quality  versus  Point  wise  cell  skewness. 
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2ET.S  CGS^USTRwirr^-i?'00  CG1  <STE  W'*h  “‘'“0  CG21  |STR  «“  ■“•92» 


Figure  17.  Lift  anti  drag  coefficients  of  CGI,  CG21,  and  CG32. 


Figure  18.  Hybrid  mesh  of  CG42  with  75  prismatic  layers  near  the  wall. 
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Figure  19.  Solution  errors  against  number  of  prism  layers  in  the  hybrid  grids. 


Figure  20,  Lift  and  drag  coefficients  of  CG42  mesh. 
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Figure  21.  OGl  grid  overview. 
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(a)  density  residuals 


£0000 


Figure  22.  Effects  of  initial  grid  spacing  (AT.)  on  the  solution  convergence  and  errors  of  open  airfoil  grids 


Figure  23-  Lift  and  drag  coefficients  of  GGi,  OG21 
hybrid  fine  grid  (L7M  cells)  of  Ref.  2i. 


and  OG32.  The  coefficients  were  compared  with  the 
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(a)  Closed  inlet 


(b)  Open  inlet 


Figure  24 , 


Effects  of  initial  grid  spacing  (^*1)  on  the  solution  errors  of  closed-  and  open-inlet  grids. 


(a)  OG5  (b)  OG12 

Figure  25*  GGI2  (with  148,262  cells)  grid  overview  compared  to  OG5  mesh  (with  410,260  cells). 


(a)  density  residuals  (h)  errors 

Figure  26.  Effects  of  cavity  cell  numbers  on  the  solution  convergence  and  errors  of  open  airfoil  grids . 
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Figure  27,  Pressure  plots  of  open-  and  closed-inlet  at  a  =  &5°. 


(a)  OG20 


(b)  OG22 


Figure  28,  Overview  of  GG2G  and  QG22  meshes. 


(b)  drag  coefficient 


Figure  29,  Lift  and  drag  coefficients  of  OG5,  OG12,  and  OG20> 
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Figure  30,  OG40  mesh  overview. 
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(a)  lift  coefficient  (b)  drag  coefficient 

Figure  31.  Lift  and  drag  coefficients  predicted  by  the  hybrid  OG40  mesh. 
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Vortical  Flow  Prediction  of  the  AVT-183  Diamond 

Wing 

Mehdi  Ghoreyshi* *  Krzysztof  Ryszka" 

Russell  M.  Cummings1  and  Andrew  J.  Lofthouse§ 

High  Performance  Computing  Research  Center  r  Ik  S.  Air  Force  Academy 
USA  F  A  ca de my *  Colo rado  S0S4  0-  64  00 


The  objective  of  this  paper  is  to  assess  the  potential  and  limitations  of  current  practice 
in  computational  fluid  dynamic  modeling  for  predicting  vortical  flowfields  over  a  generic 
53-degree  swept  diamond  wing  with  blunt  tips*  This  whig  was  designed  under  SI  O  AVI 
Task  Group  183  and  is  based  on  the  SACCON  UCAV  geometry,  which  is  a  lambda  wing 
with  complex  span-wise  distributions  of  thickness  and  leading  edge  radius  and  with  a  linear 
twist  onboard  of  the  first  trailing  edge  break*  The  SACCON  wing  trailing-edge  was  swept 
forward  by  26,5-degrees  to  form  a  diamond- shaped  planform,  This  new  wing  also  has  a 
constant  airfoil  section  of  NACA  64A-0O6  with  a  leading  edge  radius  of  0,264  in  percent 
chord  and  has  no  twist  angle,  CFD  simulations  were  performed  for  an  angle  of  attack  sweep 
at  Mach  number  of  0.15  and  Reynolds  number  of  2,7x  10°  based  on  the  mean  aerodynamic 
chord  to  match  experiments.  CFD  simulations  were  run  with  different  turbulence  models, 
that  close  the  Reynolds  averaged  Navier- Stokes  equations,  and  with  a  limited  assessment  of 
delayed  detached  eddy  simulations*  The  wind  tunnel  experiments  of  the  diamond  wing  were 
carried  out  in  the  Institute  of  Aerodynamics  and  Fluid  Mechanics  of  Technische  Universitat 
Mlinc hen,  Germany  and  include  aerodynamic  lift,  drag,  and  pitch  moment  measurements 
as  well  as  span- wise  pressure  distributions  at  different  chord- wise  locations.  This  data  set 
is  used  to  validate  CFD  results.  The  results  presented  demonstrate  that  the  pitch  moments 
predicted  by  the  SARC  turbulence  model  provide  a  better  match  to  experimental  results 
than  the  SA  model  at  moderate  angles  of  attack.  However,  at  high  angles  of  attack,  all 
pitch  moment  predictions  are  off.  The  flow  visualization  results  show  that  a  leading-edge 
vortex  is  formed  above  the  upper  surface  of  the  wing  at  an  angle  of  attack  about  eight 
degrees.  This  vortex  will  become  larger  and  stronger  as  the  angle  of  attack  is  increased. 
With  increasing  angle  of  attack,  the  vortex  starting  point  moves  upstream  and  the  vortex 
core  moves  inboard  towards  the  wing  root*  The  results  also  show  that  for  the  range  of 
angles  of  attack  considered  in  this  work*  the  flow  over  the  diamond  wmg  is  steady*  because 
the  vortices  are  attached  and  stationary  at  each  angle  of  attack. 


*  Senior  Aerospace  Engineer,  AIAA  Senior  Member 

*Cadct  of  Aeronautical  Engineering 
t Professor  of  Aeronautics,  AIAA  Associate  Fellow 
§  Director.  AIAA  Senior  Member 
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Nomenclature 


Cl 

lift  coefficient,  L/q^S 

Cm 

pitch  moment  coefficient,  Mfq. 

C 

mean  aerodynamic  chord,  m 

Cr 

chord  length  at  wing  root,  m 

L 

lift  force,  N 

M 

Mach  number,  V/a 

M 

pitch  moment,  N-m 

c 

dynamic  pressure,  Pa,  pV2  / 2 

Re 

Reynolds  number,  pVc/p 

s 

half  span,  m 

t 

time,  s 

V 

freestream  velocity,  m/s 

x,y,z 

aircraft  position  coordinates 

Greek 

a 

angle  of  attack,  rad 

8 

side-slip  angle,  rad 

/* 

air  viscosity 

I.  Introduction 

Today’s  modern  fighter  aircraft  typically  employ  a  delta  wing  configuration  to  reduce  the  wave  dra-* 
known701ih  S,PefS'  f  aer°dynamiC  performances  of  these  delta  wings  are  vastly  different  from  those 

the  attache  ‘fl  7 *  ?  ™gS  "  f”  latter  case-  the  lift  Besses  linearly  with  angle  of  attack  in 

,  *  d  fl,°  .  [egl0n  a!ld  '-hen  sharply  decreases  in  the  post-stall  region,1  The  lift  of  a  sharp  delta  wing 

to  tit  ZZi  TtT’  T  ^  fn  fDgle  °f  attack  °f  just  few  ^ere  is  an  additional  lift  force 

7.  |  W  flow  hft  Which  makes  the  lift-curve  slope  nonlinear.  This  additional  lift,  often  called  the 

77  f  i  V7iC  2  7WS  f°rmed  ab°Ve  the  Wlng‘2  The  flow  separation  point  of  a  sharp  delta 
'  .  '  .  e  eadnig  edge.  This  creates  a  strong  shear  layer  along  the  wing  edges  which  will  then 

up  into  a  pair  of  counter-rotating  vortices  over  the  upper  surface  on  the  two  halves  of  the  wing  3  These 
vortices  are  called  leading-edge  vortices  and  their  structures  depend  on  the  wing’s  sweep  angle,  leading-edge 
geometry,  wmg  thickness,  and  freestream  conditions,4 

ring’ the Sapara*ion  POintS  are  fixed  at  the  ]eadill§  edSe for  a  considerable  range 
°f  7  7  i  f  ,'  r  r  StrellgtllS  aild  VOrtex  lift  tIierefore  wil1  increase  with  the  increase  in  the  angle 
of  attack.  At  higher  angles,  however,  vortices  experience  an  abrupt  transformation  called  the  vortex  break- 

771  which  is  an  asymmetric  phenomenon.  In  a  vortex  breakdown,  the  axial  velocitv  component  suddenly 
decelerates  and  the  swirl  component  of  the  mean  velocity  decreases  due  to  the  vortex  core  expansion  0  The 
asymmetric  vortex  breakdown  conditions  result  in  additional  moments  in  pitch,  yaw.  and  roll  with  magni¬ 
tudes  as  large  or  even  larger  than  those  obtained  from  traditional  control  surface  deflections  The  vortex 

cnnS?i°^r7y  T  ?  ?n  C‘rSeS  “  Pitching  moment’  loss  of  i,ft-  a!ld  buffeting.7  The  stability  and 
control  (S&C)  and  structural  analysis  of  delta  wings  are  therefore  highly  dependent  on  the  ability  to  observe 
and  accurately  calculate  the  principle  characteristics  of  vortical  flows. 

While  the  vortical  A™  behavior  over  slender  wings  with  sharp  leading  edges  has  been  studied  extensively 
lor  over  last  few  decades,  ■  •  much  less  is  known  about  the  formation  of  vortices  over  wings  with  lower 
sweep  angles  and  blunt  leading  edges.  These  type  of  wings  are  often  incorporated  in  the  designs  of  unmanned 

Zentrf  7  ,  7CAT  ,F°r  tlMSe  Wi“gS’  th*  ^  flow  structure  is  very  implicated  and 
depends  heavily  on  the  leading  edge  bluntness  and  the  wing  sweep  angle.  For  delta  wings  with  blunt  tips 

the  leading  edge  separation  point  is  also  very  sensitive  to  the  boundary  layer  changes.12  The  blunt  tip 

777  77  *77  TIT  fffects  on  the  stab0it.V  aild  c01'1*™1  characteristics  of  a  maneuvering  aircraft 
at  high  angles  of  attack  The  objective  of  this  work  is  to  assess  the  potential  and  limitations  of  current 
practice  in  computational  fluid  dynamics  (CFD)  modeling  for  predicting  vortical  flowfields  over  a  generic 
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5 3- degree  swept  diamond  wing  with  blunt  tips. 

The  diamond  wing  considered  in  this  paper  was  designed  under  STO  AVT  Task  Group  183.  I  he  wing 
plan  form  is  based  on  the  SACCON  UCAV  geometry.  The  SACCON  wing  trailing-edge  was  swept  forward  by 
26.5-degrees  to  form  a  diamond- shaped  planfonn.  The  new  wing  has  a  constant  airfoil  section  of  NACA  64 A- 
006  with  a  leading  edge  radius  of  0.264  in  percent  chord.  The  wind  tunnel  experiments  of  the  diamond  wing 
were  carried  out  in  the  Institute  of  Aerodynamics  and  Fluid  Mechanics  of  Technische  Uni  vers  it  at  Miinchen, 
Germany  and  include  aerodynamic  lift,  drag,  and  pitch  moment  measurements  as  well  as  span- wise  pressure 
distributions  at  different,  chord-wise  locations.  This  data  set  is  used  to  validate  CFD  results. 

This  work  is  organized  into  four  major  categories.  First  a  survey  of  the  literature  and  theories  pertaining 
to  the  delta  wing  vortical  flows  is  provided.  The  flow  solver  is  then  briefly  described.  Next,  the  test  case  is 
presented  and  the  experimental  setup  is  detailed.  Finally,  the  simulation  results  will  be  discussed. 

II.  Vortical  Flows 

The  vortical  flow  behavior  over  delta  wings  is  extremely  complex  and  differs  substantially  from  sharp  to 
blunt  tips  and  from  nonslender  to  slender  wings.  Over  the  past  few  decades,  most  research  was  done  on 
the  vortices  over  slender  sharp-edged  delta  wings.13  Some  early  works  in  this  field  were  reviewed  by  Sun 
and  most  recently  by  Mitchell  et  al.4  In  one  of  the  earliest  studies,  for  example,  Jensen  (1948)lj  studied 
the  low-speed  flowfields  and  the  lift  and  moment  characteristics  of  a  sweptback  wing  and  a.  delta  wing,  both 
with  a  65°  sweptback  leading  edge  and  a  double  wedge  symmetrical  airfoil  section.  His  results  showed  that 
two  strong  vortices  are  formed  over  the  upper  surfaces  of  both  wings  for  angles  of  attack  as  low  as  10  .  A 
sketch  of  these  vortices  is  shown  in  Figure  1  for  the  delta  wing  case.  Jensen's  work  also  showed  that  tire  lift 
curve  slope  increases  for  angles  of  attack  above  ten  degrees. 

It  is  well  known  that  a  sharp  leading  edge  causes  the  boundary  layer  to  separate  at  t  he  leading  edge 
resulting  a  free  shear  layer.  For  a  slender  wing,  the  separated  shear  layers  roll  up  into  a  pair  of  counter¬ 
rotating  vortices  over  the  upper  surface  on  the  two  halves  of  the  wing.  The  shear  layer  may  exhibit  instabilities 
that  increase  the  vortical  substructures  and,  therefore,  the  primary  vortex  increases  in  both  size  and  strength 
as  it  extends  downstream.  Ornberg16  in  1954  showed  that  these  leading  edge  vortices  are  stationary  and 
form  a  low  pressure  region  over  the  upper  surface  that  will  increase  the  lift.  Notice  that  this  additional 
lift,  comes  at  the  expense  of  a  drag  penalty  due  to  loss  of  leading-edge  suction,17  The  leading  edge  vortices 
allow  the  onset  of  stall  to  be  delayed  to  higher  angles  as  well.  For  steady  and  in  vise  id  flow,  the  vortex  lift 
can  be  approximated  to  some  extent  by  the  leading  edge  suction  analogy,18  linear  slender  wing  theory, 
or  detached  flow  methods.20  The  Polhamus's  leading  edge  suction  analogy  is  probably  the  most  widely 
applicable  method  to  estimate  the  vortex  lift  of  different  planfonns  at  different  flight  speeds.  This  analogy 
assumes  that  the  vortex  lift  has  the  same  magnitude  as  the  potential-flow  force  which  is  lost  because  of  the 
separation  at  the  sharp  leading  edge.21  The  method  lias  provided  good  results  ior  the  attached  vortices. ‘s 

As  the  angle  of  attack  increases,  the  core  of  leading  edge  vortices  will  move  inboard  and  secondary 
vortices  can  be  formed  below  the  main  vortices.22  Ornberg16  in  1954  and  Marsden  and  his  colleagues21  in 
1958  reported  the  existence  of  these  secondary  vortices.  They  found  that  the  secondary  vortices  form  below 
and  outboard  the  main  wing  vortices  with  opposite  circulation.  The  leading-edge  vortices  can  separate  near 
the  surface  because  of  the  adverse  pressure  gradients  in  the  spanwise  direction.23  These  separated  flows 
may  then  roll  up  to  form  smaller  secondary  vortices  of  the  opposite  sign  of  vorticity  which  tend  to  move  the 
primary  vortices  inboard  and  upward.' 

Leibovich24  showed  that  in  a  slender  delta  wing,  each  leading-edge  vortex  core  is  approximately  one  core 
diameter  above  the  lifting  surface  of  the  wing.  Over  a  nonslender  wing,  however,  the  vortices  are  formed 
much  closer  to  the  surface  of  the  wing,  therefore,  the  interaction  of  the  secondary  vortices  with  the  shear 
layer  of  the  primary  vortex  causes  the  primary  vortex  to  split  into  a  dual  primary  vortex  structure.  "J  Since 
the  cores  of  vortices  over  nonslender  wing  are  much  closer  to  the  surface  and  its  boundary  layer,  the  flow 
over  nonslender  wings  is  more  sensitive  to  Reynolds  number  than  flows  over  slender  wings.* 1  The  vortical 
flow  separation  also  occurs  at  lower  angles  of  attack  compared  with  a  slender  delta  wing. 

The  leading-edge  vortices  have  significant  effects  on  the  aerodynamic  behavior  of  aircraft;  they  form 
regions  of  high  vorticity  and  low-pressure  over  the  upper  wing  surface  that  cause  a  nonlinear  increase  in 
lift  until  the  maximum  attainable  lift  with  attached  flow.8  At  high  angles  of  attack,  the  vortex  structures 
change  dramatically,  leading  to  vortex  breakdown,  which  is  known  to  cause  nonlinear  aerodynamic  behavior. 
Vortex  breakdown  is  characterized  by  a  sudden  deceleration  of  the  axial  flow*  in  the  vortex  core,  a  decrease  in 
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circumferential  velocity,  and  an  increase  in  the  vortex  size  and  turbulent  dissipation.26  The  onset  of  vortex 
breakdown  starts  from  the  trailing  edge  of  the  wing  and  moves  forward  on  the  wing  with  increasing  aiwle 
of  attack  until,  at  sufficiently  high  angle  of  attack,  the  flow  is  dominated  by'  the  periodic  wake  shedding 
encountered  over  bluff  bodies.2'  Note  that  a  slender  wing  attains  its  maximum  lift  and  stalls  approximately 
when  the  onset  of  vortex  breakdown  crosses  the  trailing  edge,  while  wings  with  smaller  sweep  angles  stall 
when  the  onset  of  vortex  breakdown  moves  towards  the  wing  apex.26 

One  of  the  methods  to  reduce  the  vortex  drag  is  using  a  rounded  leading  edge.28  While  the  leading  edge 
suction  is  (normally)  lost  in  the  sharp  edges,  a  large  portion  of  the  leading-edge  suction  is  restored  in  a  round 
e  ge  an  ence  tin’  drag  force  becomes  smaller.  However,  the  blunt  tips  fundamentally  change  the  vortical 
flow  behavror  and  structure.  On  delta  wings  with  rounded  leading  edges,  the  flow  separation  point  is  not 
fixed  at  the  leading  edge  and  depends  on  the  Reynolds  number  and  the  leading  edge  curvature.  The  origin 
o  t  re  leading  edge  vortex  is  at  the  delta  wing  apex  for  a  sharp-edged  wing  but  it  is  further  downstream  for 
a  rounc  e  tip.  lie  leading-edge  vortices  from  sharp  edges  are  steady  for  a  wide  range  of  angles  of  attack. 

n  y  if  the  angle  of  attack  increases  greatly,  these  vortices  become  unsteady  and  separate  starting  from  the 
trailing  edge  of  the  wmg.  The  flow  at  blunt  nose,  however,  separates  at  moderate  to  high  angles  of  attack 
starting  from  an  outboard  and  aft  location  near  the  trailing  edge29 

III.  CF D  Solver 

Cobalt  solves  the  unsteady,  three-dimensional,  compressible  Navier-Stokes  equations  in  an  inertial  ref¬ 
erence  frame  .Arbitrary  cell  types  in  two  or  three  dimensions  may  be  used;  a  single  grid  therefore  can  he 
composed  of  different  cell  types.30  In  Cobalt,  the  Navier-Stokes  equations  are  discretized  on  arbitrary  grid 
topologies  using  a  cell-centered  finite  volume  method.  Second-order  accuracy  in  space  is  achieved  using  the 
exact  Riemann  solver  of  Gottlieb  and  Groth  31  and  least  squares  gradient  calculations  using  QR  factoriza¬ 
tion.  To  accelerate  the  solution  of  the  discretized  system,  a  point-implicit  method  using  analytic  first-order 
mviscid  and  viscous  Jacobians  is  used.  A  Newtonian  sul> iteration  method  is  used  to  improve  the  time  ac¬ 
curacy  of The  point-implicit  method.  Tomaro  et  al.32  converted  the  code  from  explicit  to  implicit,  enabling 
Cou rant- Ftie d richs- Le wy  (CFL)  numbers  as  high  as  10°.  In  Cobalt,  the  computational  grid  can  be  divided 
m  o  group  of  cells,  or  zones,  for  parallel  processing,  where  high  performance  and  scalability  can  be  achieved 
°n  ^  34  oUSaindS  °fprocessors; '  Some  available  turbulence  models  in  Cobalt  are  the  Spalart-AUmaras 
SJ.K,  Spaiart-AOmaras  with  Rotation  Correction  (SARC),35  and  Delayed  Detached-eddy  simulation 
[DDhjb)  with  SARC. 

IV.  Test  Case 

The  diamond  wing  used  for  this  work  was  designed  within  the  NATO  RTO  task  group  AVT-183  The 
objective  of  this  task  group  was  reliable  prediction  of  separated  flow  onset  and  progression  for  air  and  sea 
vehicles.  The  AV  T-183  wing  is  based  on  the  SACCON  UCAV  geometry  of  AVT-161  research  program37  but 
it  retains  a  simple  and  less  complicated  design  compared  to  SACCON.  In  more  details,  the  SACCON  has 
a  lambda  wmg  planform  with  a  leading  edge  sweep  angle  of  53°  as  shown  in  Fig.  2.  The  main  sections  of 
the  SACCON  model  are  the  fuselage,  the  wing  section,  and  wing  tip.  The  configuration  is  defined  by  three 
different  profiles  at  the  root  section  of  the  fuselage,  two  sections  with  the  same  profile  at  the  inner  wing 
lor ming  the  transition  from  the  fuselage  to  wing  and  the  outer  wing  section. 

1  he  spanwise  distributions  of  thickness  in  percent  chord  and  leading-edge  radius  of  SACCON  are  shown 
m  Figure  3.  According  to  this  figure,  the  leading  edge  radius  is  sharp  at  the  root  chord  and  then  increases 
in  the  spanwise  direction  up  to  the  intersection  between  fuselage  and  wing  and  then  decreases.  Figure  3 
also  shows  that  the  thickness  ratio  in  spanwise  direction  decreases.  Finally,  the  outer  wing  section  profile  is 

twisted  by  5°  around  the  leading  edge  to  reduce  the  aerodynamic  loads  and  shift  the  onset  of  flow  separation 
to  higher  angles  of  attack. 

oJ^rbmati°n  ,°fEW?P  angle  With  sharP-bIunt  leading  edges  and  the  twist  angle  makes  the  flow  around 
SACCOIS  very  complicated  and  difficult  to  predict  even  for  static  conditions.38-39-40  This  is  due  to  formation 
a  an  apex  vortex  and  a  tip  vortex  which  are  formed  by  the  onset  and  progression  of  the  flow  at  the  sharp  and 
:>lunt  leading  edges  at  moderate  angles  of  attack.  At  higher  angles  of  attack  these  vortices  become  stronger 
and  the  onset  of  the  tip  vortex  moves  inboards.  Over  a  small  range  of  angle-of-attack,  two  vortices  merge 
and  coalesce  to  form  much  a  stronger  and  deeper  vortex  which  suddenly  changes  the  pitch  moment  values. 
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At  higher  angles,  vortices  breakdown  and  again  cause  severe  pitch  moment  changes. 

For  all  above  reasons,  although  CFD  results  of  AVT-161  match  each  other  at  small  angles  of  attack  but 
they  predicted  different  values  at  moderate  to  high  angles  of  attack;  besides  none  ol  codes  could  match  the 
SACCON  experimental  data  for  the  full  range  of  angles  of  attack.  To  follow  up  the  vortex  flow  studies 
of  AVT-161  and  gain  information  to  understand  why  some  CFD  codes  miss  or  match  the  leading-edge 
separation  aerodynamics  of  interest,,  the  work  of  AVT-183  has  been  focused  on  a  new  geometry  with  simpler 
design  complexity  than  the  SACCON.  The  new  wing  has  a  constant  radius  blunt  leading  edge  m  spanwise 
direction.  The  thickness  ratio  is  also  constant  and  the  wing  is  non-twisted.  In  addition,  the  trailing  edge 
has  swept  forward  26.5-degrees  to  form  a  diamond  wing.  This  design  has  been  compared  with  the  SACCON 
geometry  in  Fig.  4. 

Figure  5  shows  the  diamond  wing  layout  used  in  the  wind  tunnel  experiments.  The  wind  tunnel  model 
was  places  on  a  peniche  with  height  of  9  cm  to  minimize  the  boundary  layer  effects  from  the  wing  tunnel 
wall.  Figure  5  shows  that  the  diamond  wing  wind  t  unnel  model  has  a  root  chord  of  1.2  m,  slightly  bigger 
than  SACCON.  The  wing  half  span  is  approximately  0.65m.  The  moment  reference  point  (aWp  is  located  at 
0.49  m  from  apex.  A  constant  NACA64A006  airfoil  section  was  used.  The  airfoil  shapes  at  some  spanwise 

locations  are  shown  in  Figure  6.  , 

A  hybrid  RANS  mesh  was  generated  for  this  model.  The  mesh  was  generated  in  two  steps.  In  the  first 
step,  the  mviscid  tetrahedral  mesh  was  generated  from  a  clean  configuration  using  the  ICEM-CFD  code. 
The  inviscid  mesh  was  then  used  as  a  background  mesh  by  TRITET41,42  which  builds  prism  layers  using  a 
frontal  technique.  TRITET  rebuilds  the  inviscid  mesh  while  respecting  the  size  of  the  original  inviscid  mesh 
from  ICEM-CFD.  Mesh  overview  is  shown  in  Figure  7.  The  grid  is  symmetric  configuration  and  contains 
5.4  million  points  and  13.2  million  cells. 


V.  Wind  Tunnel  Experiments 

A  half  model  of  the  diamond  wing  was  built  for  wind  tunnel  testing.  The  model  layout  is  shown  in 
Figure  5.  The  experiments  of  the  diamond  wing  were  carried  out  at  the  Technische  Universitat  Miinchen 
subsonic  wind  tunnel.  Figures  8  shows  the  wing  model  with  a  peniche  in  the  test  section  of  this  wind  tunnel. 
The  wind  tunnel  is  operable  in  a  closed-  and  open-circuit  mode  with  a  suction  configuration.  lest,  section 
area  is  l.SOmx 2.40m  with  a  length  of  4.8  meters.  The  wind  tunnel  maximum  speed  is  about  65  m/sec. 

All  experiments  were  conducted  for  an  angle  ol  attack  sweep  at  Mach  number  ol  0,15  and  Reynolds 
number  of  2.7x  10®  based  on  the  mean  aerodynamic  chord.  An  external  six-component,  balance  was  used 
to  measure  the  forces  and  moments  acting  on  the  model  only  and  not  the  peniche.  A  stereo  pai  tide  image 
velocimetry  (PIV)  technique  was  also  used  to  measure  flow  field  properties  and  detect  vortical  flows.  The 
PIV  technique  was  done  with  a  laser  and  two  cameras.  The  camera  and  laser  positions  aic  highlighted 
in  Figure  8.  For  each  test  with  different  camera  angles,  about  400  images  were  captured.  These  images 
were  processed  using  LaVisiom  DaVis  software.  The  experimental  data  include  averaged  lift,  drag,  and  pitch 
moment  coefficients  for  all  test  runs  over  the  range  of  tested  angles  of  attack.  Besides,  the  surface  pressure 
coefficients  were  measured  for  a  number  ol  slices  at,  different  chord-wise  locations.  These  slices  are  shown  in 
Figure  9. 


VI.  Results  and  Discussion 

Steady  and  unsteady  simulations  of  the  diamond  wing  were  carried  out  with  different  turbulence  models. 
All  simulations  retain  second-order  spatial  accuracy  and  use  an  implicit  temporal  operator  to  advance  the 
flow  solution  in  time  However,  the  steady-state  cases  are  first-order  accurate  in  time  and  were  computed 
with  a  CFL  of  one-million.  Unsteady  computations  were  performed  using  second-order  temporal  accuracy 
and  three  Newton  sub-iterations.  The  effects  of  the  time  step  on  the  predictions  are  evaluated  for  unsteady 

simulations.  t 

Turbulence  models  considered  include  the  Sp  al  art  -  A 1 1  ni  aras  (SA)  model,  Spalar  t-  A  llmaras  with  Rota- 
lion  Correction  (SARC),  and  Delayed  Detached-eddy  simulation  (DDES)  with  SARC.36  Note  that  the 
DDES  formulation  will  improve  predictions  of  flows  with  massive  separation  which  are  time-dependent  flows. 
Therefore,  all  DDES  simulations  were  performed  using  second-order  temporal  accuracy  and  three  Newton 
sub- iterations. 
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I  he  computed  overall  lift,  drag,  and  pitch  moment  coefficients  using  the  SA  turbulence  model  Eire  shown 
in  igure  10  and  compared  to  experimental  data.  These  simulations  are  steady-state  type  and  correspond 
to  an  angle-of-attack  sweep  from  zero  to  20=  with  a  0.5  degrees  increment  for  angles  between  8  to  16  degrees 
and  one  degree  increment  for  all  other  angles.  These  angles  are  selected  allowing  us  to  compare  CFD 
predictions  with  static  pressure  measurements  in  the  wind  tunnel  at  each  of  these  angles.  Two  CFD  data 
sets  are  presented  in  Figure  10;  the  first  dataset  are  simulated  for  1,500  iterations  for  all  angles  of  attack;  the 
second  dataset  corresponds  to  simulations  that  run  for  more  iterations  (3,000)  at  some  high  angles  of  attack. 
Figure  10  shows  that  the  computed  overall  lift,  drag,  and  pitch  moment  coefficients  from  both  simulations 
ma.  cli  very  well,  indicating  that  the  force  and  moment  coefficients  have  already  reached  their  final  values  at 
iteration  1,500  and  they  do  not  change  further  in  time. 

Figure  10  shows  that  predictions  of  CFD  with  the  SA  turbulence  model  match  very  well  with  experiments 

VT  cS  °  attaCk  but  they  d°  “0t  match  at  higher  ailgles’  excePt  the  drag  force.  At  high  angles  of 
a'dC/  ,  c‘ata  overestimate  and  pitch  moment  predictions  underestimate  the  measurements;  also  the 
experiments  show  that  there  is  a  break  in  the  pitch  moment  curve  around  8°  angle  of  attack;  CFD  using  the 
SA  model,  however,  shows  the  break  in  pitch  moment  curve  at  a  much  higher  angle  (around  12°). 

Results  of  steady-state  simulations  with  the  SARC  turbulence  models  are  compared  witli  experiments 
and  predictions  obtained  from  the  SA  model  in  Figure  11.  These  simulations  were  first-order  accurate  in 
irne  and  were  computed  with  a  CPL  of  one-million;  however  t  hey  ran  for  6,000  iterations  to  allow  forces  and 
moments  reaching  t  heir  final  values.  Figure  11  shows  that  the  SARC  and  SA  turbulence  model  predictions 
are  in  good  agreement  with  each  other  and  experiments  at,  small  angles  of  attack.  However,  at  moderate  to 
ug  i  angles  predictions  from  these  turbulence  models  are  different;  the  SARC  drag  values  are  slightly  larger 
than  the  SA  model  as  well  as  experiments;  the  lift  coefficients  from  the  SARC  model  are  slightly  smaller 
than  predicted  lift  from  the  SA  model  and  become  closer  to  the  experiments  for  angles  of  attack  above  15 
egrees^  1  he  SARC  pitch  moment  values  agree  well  with  experiments  at  moderate  angles  of  attack  but  they 
are  o  for  angles  above  15  degrees.  The  SARC  model  accurately  predicts  the  pitch  moment  break  seen  in 
the  experiments  around  8°  angle  of  attack. 

i  FilgUJ!  o2r,ShOWf.  Ule  <,ffectS  °f  the  time  steP  size  i!1  the  unsteady  SARC  computations  compared  with  the 
steady  SARC  predictions.  Unsteady  computations  were  performed  using  second-order  temporal  accuracy 
and  three  Newton  sub-iteration  with  time  step  sizes  of  2e-4  and  4e-5  seconds.  These  unsteady  simulations 
were  ran  for  6,000  iterations  and  the  overall  forces  and  moments  were  averaged  for  the  last  3,000  iterations, 
rhe  resuits  in  Figure  12  reveal  that  discrepancies  between  the  steady  and  unsteady  predictions  are  verv 
small,  this  indicates  that  the  flow  is  steady  and  attached  over  the  surfaces  in  this  angle  of  attack  range. 

•  lgrSo  l3  and  £«SeUt, the  fl°W  visualization  the  diamond  wing  at  different  angles  of  attack  bv 

using  the  SA  and  SARC  turbulence  models,  respectively.  Visualization  of  vortices  over  the  upper  surfaces 
was  tone  .y  using  iso-surfaces  of  vorticity;  the  surface  and  vortices  in  these  figures  were  colored  with  the 
pressure  coefficient  ranging  from  -3  to  2.  These  figures  show  that  visualization  plots  for  a  =4°  are  very 
similar  for  the  SA  and  SARC  models.  At  this  angle,  the  flow  separates  at  the  wing-tip  near  trailing  edge 
due  to  considerable  pressure  difference  between  the  lower  and  upper  surfaces.  This  separated  flow  then  rolls 
up  into  a,  small-size  tip  vortex  which  can  be  seen  in  Figures  13(a)  and  14(a).  These  figures  show  that  at  the 
wing  leading  edge,  there  is  a  positive  spanwise  pressure  gradient  towards  the  center  of  wing  The  pressure 
gradient  becomes  larger  at  the  wing  tip. 

At  six  degrees  angle  of  attack,  the  tip  vortex  becomes  stronger  as  shown  in  Figures  13(b)  and  14(b) 
However,  the  size  of  the  vortex  predicted  by  the  SARC  turbulence  model  is  larger  than  the  one  predicted  by 
,  ie  .  .  model.  1  lie  vortex  predicted  from  SARC  is  also  located  further  upstream.  Notice  that  the  pressures 
at  the  leading  edge  are  smaller  than  four  angle  of  attack  case;  the  spanwise  pressure  gradient  is  therefore 
larger  at  six  degrees  angle  of  attack.  Figures  13  and  14  show  that  at  about  8=  angle  of  attack,  a  leading  edge 
voitex  is  formed  above  the  wing*  This  leading-edge  vortex  originates  from  about  0,4  and  0.6  cr  (the  chord 
length  at  the  wing  root)  for  the  SARC  and  SA  models,  respectively.  Again  the  SARC  predicted  vortex  is 
bigger  and  stronger  than  the  vortex  predicted  by  the  SA  model.  These  leading-edge  vortices  merge  with  the 

tip  vortex.  A  low  pressure  region  is  formed  above  the  upper  surface  around  the  vortex  core.  Vortices  also 
increase  in  size  as  they  extend  downstream. 

As  angle  of  attack  is  further  increased,  the  leading-edge  vortices  become  stronger  and  bigger;  the  vortex 
s tarring  point  also  moves  upstream.  For  example  at  ten  degrees  angle  of  attack,  the  starting  points  of  vortices 
predicted  by  the  SA  and  SARC  models  are  at  about  0.5  and  0.3  c,,  respectively.  The  center  of  the  vortex 
also  moves  inboard  with  increasing  angle  of  attack. 
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More  details  of  the  flows  around  the  diamond  wing  can  be  found  in  Figures  15-22  which  compare  pre¬ 
dicted  pressure  coefficients  by  the  SA  and  SARC  model  against  experiments  for  several  angles  of  attack. 
In  these  figures,  the  corresponding  pressure  distributions  on  the  upper  and  lower  surfaces  of  the  wing  are 
shown.  Figure  15  shows  that  at  all  chord-wise  locations,  considered  in  this  work,  CFD  predictions  from  both 
turbulence  models  match  well  with  each  other  and  experiments.  In  all  plots  of  this  figure,  moving  towards 
the  wing  tip  (y/s  =  1)  a  positive  and  a  negative  pressure  gradient  exist  at  the  lower  and  upper  surface, 
respectively.  The  pressure  curves  for  a  =4°  are  smooth  with  a  minimum  pressure  at  the  wing  tip.  figure  15 
shows  that  no  vortex-type  flow  is  present  for  the  range  of  chordwise  locations  considered. 

Figure  16  compares  the  pressure  distribution  plots  at  a  =6°.  The  pressure  differences  between  the  upper 
and  lower  surfaces  are  larger  than  plots  at  a  =4°.  CFD  pressure  values  match  with  experiments  almost 
everywhere.  The  largest  discrepancies  in  plots  can  be  seen  at  0.6cr.  At  this  location,  the  SARC  model 
predicts  a  negative  pressure  region  in  the  upper  surface  near  the  wing  tip  due  to  the  tip  vortex  presence;  the 
vortex  core  is  located  about  y/s  =  0.93.  However,  SA  and  experiment  results  do  not  predict  the  existence  of 
such  a  vortex  at  this  location  and  angle.  Figure  16  shows  that  at  0.6cr,  there  is  a  positive  span  wise  pressure 
gradient  towards  the  wing  tip  at  the  vortex  boundaries.  The  lower  surface  pressure  distributions  at  this 
angle  of  attack  are  smooth  and  have  no  interesting  flow  features  present. 

The  pressure  differences  between  the  upper  and  lower  surfaces  at  a  =8°  even  become  larger.  The  lower 
surface  pressure  distributions  show  no  interesting  flow  feature  again  as  shown  in  Figure  17.  However,  the 
upper  surface  lias  a  strong  negative  dip  at  some  chordwise  locations  due  to  the  presence  of  a  leading  edge 
vortex.  The  SARC  model  exhibits  the  vortex  presence  even  at  0.395cr.  As  flow  moves  downstream,  the 
predicted  vortex  becomes  bigger  (the  dip  becomes  wider),  however  the  pressure  values  become  less  negative. 
The  SA  predicts  a  vortex  above  the  upper  surface  as  well.  This  vortex  is  seen  only  at  0.6cr.  Experiment 
results  also  exhibit  a  vortex  starting  around  0.5cr.  Figure  17  shows  that  pressure  distributions  predicted  by 
CFD  do  not  exactly  match  with  experiments  around  vortices. 

As  angle  of  attack  is  increased  to  a  =10°,  the  SARC  predicted  vortex  is  moved  further  upstream,  at 
about  0.295cr  as  shown  in  Figure  18.  The  vortex  core  is  also  moved  inboard.  For  example,  Figure  17 
show's  that  for  a  =8°  at  0.395cr,  the  vortex  core  is  located  about  y/s  =  0.95:  the  vortex  core  is  located  at 
about  y/s  =  0.85  for  a  =10°.  Experiment  results  exhibit  that,  a  vortex  is  formed  starting  at  0.395cr.  The 
SA  model,  however,  predicts  the  vortex  starting  point  at  0.5cr.  CFD  data  and  experimental  results  show' 
that  the  predicted  vortex  becomes  larger  in  size  and  the  pressure  values  become  less  negative  as  the  flow 
moves  upstream.  Also,  a  positive  spanwise  pressure  gradient  is  predicted  towards  the  wing  tip  at  the  vortex 
boundaries  in  the  simulation  and  experimental  results. 

At  12  degrees  angle  of  attack,  experimental  results  shows  a  vortex  starting  at  0,395cr;  the  SARC  and  SA 
vortices  now  originate  from  0.295cr  and  0.395cr,  respectively.  The  vortices  become  bigger  and  induce  more 
negative  pressure  regions  as  shown  in  Figure  19.  For  example,  at  0.6cr  the  vortex  diameter  is  about  0.4  ol 
half-span  length.  As  the  angle  of  attack  is  further  increased,  vortices  predicted  from  CFD  and  experiments 
become  much  bigger.  At  high  angles  or  attack,  the  pressure  distribution  becomes  more  uniform  across  the 
vortex  at  upstream  locations. 

Final  results  compare  CFD  pedictions  using  DDES-SARC  turbulence  model  with  experiments.  DDES 
simulations  ran  for  12,000  iterations  using  a  second-order  time  accuracy  and  a  time  step  of  4e-5  second. 
Lift,  drag,  and  pitch  moment  coefficients  were  averaged  for  the  last  6,000  iterations,  Figure  23  compares 
the  DDES-SARC  predictions  with  experiments  as  well  as  the  SARC  predictions.  The  results  show'  that  even 
predictions  from  the  DDES-SARC  model  do  not  match  with  experiments;  prediction  curves  with  the  angle 
of  attack  are  not  smooth  as  well.  Notice  that  the  flow  is  steady-state  at  these  angles  and  a  DDES  model 
probably  do  not  much  help  to  improve  predictions. 

VIT.  Conclusions 

While  the  vortical  flow  behavior  over  slender  wings  with  sharp  leading  edges  has  been  studied  extensively 
for  over  last  few  decades,  much  less  is  known  about  the  formation  of  vortices  over  wings  with  lower  sweep 
angles  and  blunt  leading  edges.  In  the  latter,  the  vortex  flow  structure  is  very  complicated  and  depends 
heavily  on  the  leading  edge  bluntness  and  t  he  wing  sweep  angle.  These  vortices  can  have  significant  effects 
on  the  aircraft  stability  and  control  characteristics. 

CFD  potential  and  limitations  for  predicting  vortical  flowfields  were  investigated  over  a  generic  53-degree 
swept  diamond  wing  with  blunt  tips.  The  diamond  wing  considered  in  this  paper  was  designed  under  STO 
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AVI'  Task  Group  183.  The  wing  platform  is  based  on  the  SACCON  UCAV  geometry.  The  SACCON  whig 
trai ling-edge  was  swept  forward  by  20. 5- degrees  to  form  a  diamond* shaped  planform.  The  new  wing  lias 
a  constant  airfoil  section  of  NACA  64A-00G  with  a  leading  edge  radius  of  0.264  in  percent  chord.  The 
vrind  tunnel  experiments  of  the  diamond  wing  were  carried  out  in  the  Institute  of  Aerodynamics  and  Fluid 
Mechanics  of  Technische  Universitat  Miinchen,  Germany  and  include  aerodynamic  lift,  drag,  and  pitch 
moment  measurements  as  well  as  span-wise  pressure  distributions  at  different  chord-wise  locations.  This 
data  set  is  used  to  validate  CFD  results. 

Turbulence  models  considered  included  the  Spalart-Allmaras  (SA)  model,  Spalart-Allmaras  with  Rotation 
Correction  (SARC),  and  Delayed  Detached-eddy  simulation  (DDES)  with  SARC.36  The  results  presented 
demonstrated  that  the  pitch  moments  predicted  by  the  SARC  turbulence  model  provided  a  better  match  to 
experimental  results  than  the  SA  model  at  moderate  angles  of  attack.  However,  at  high  angles  of  attack, 
all  pitch  moment  predictions  were  off.  The  flow  visualization  results  showed  that  a  leading-edge  vortex  is 
formed  above  the  upper  surface  of  the  wing  at  an  angle  of  attack  about  eight  degrees.  This  vortex  became 
larger  and  stronger  as  the  angle  of  attack  was  increased.  With  increasing  angle  of  attack,  the  vortex  starting 
point  moved  upstream  and  the  vortex  core  moved  inboard  towards  the  wing  root.  The  results  also  showed 
that  foi  the  range  of  angles  of  attack  considered  in  this  work,  the  flow  over  the  diamond  wing  was  steady, 
because  the  voitices  were  attached  and  stationary  at  each  angle  of  attack.  For  this  reason,  predictions  from 
the  DDES- SARC  model  also  did  not  match  with  experiments. 
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Figure  1.  Jensen’s  sketch  of  leading  edge  vortices  over  a  delta  wing  at  a  =  2Q-,15 


Figure  2.  The  S  ACC  ON  geometry. 
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Figure  3,  SACCON  thickness  and  leading-edge  radius.  This  picture  is  adapted  from  Ref  4$. 
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Figure  4.  Diamond  wing  concept.  This  picture  is  adapted  from  Ref  43, 


Figure  5.  Diamond  wing  geometry 
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Figure  6,  NACA64A006  airfoil  section  used  in  the  diamond  wing. 


Figure  7.  Diamond  wing  mesh. 
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Figure  8.  The  AVT-183  diamond  wing  at  the  Techtiische  Universitat  Munchen  wind  tunnel.44 


Figure  9.  Pressure  tap  locations* 
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(a)  Lift  coefficient 


(b)  Drag  coefficient 


Figure  10.  CFD  predictions  of  Cobalt  with  SA  turbulence  model-  Mach  number  is  0,15  with  Rc  =  2,7  10°. 
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(a)  Lift  coefficient  (b)  Drag  coefficient 


(c)  Pitch-moment  coefficient 


Figure  IX.  SARC  versus  SA  data.  In  both,  first  order  accuracy  in  time  with  one  Newton  sub-iteration  were 
used . 
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(a)  Lift  coefficient  (b)  Drag  coefficient 


(c)  Pitch-moment  coefficient 


Figure  12.  Unsteady  simulations  using  SARC  turbulence  model.  In  all  unsteady  simulations,  second  order 
accuracy  in  time  with  three  Newton  sub-iteration  were  used.  All  unsteady  coefficients  are  averaged  values  for 
last  3,000  iterations. 
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(d)  a  =  10° 


Figure  13,  Vortical  flows  over  the  AVT-183  diamond  wing  using  SA  turbulence  model. 
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(b)  Q  -  6° 


Figure  14.  Vortical  flows  over  the  AVT-183  diamond  wing  using  SARC  turbulence  model. 
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Figure  15,  CFD  pressure  tap  data  at  o  =  4°. 


Figure  16,  CFD  pressure  tap  data  at  a  =  6°. 
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x/cr= 0.1 


Figure  17,  CFD  pressure  tap  data  at  a  = 


Figure  18,  CFD  pressure  tap  data  at  a  =  10°, 
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Figure  19*  CFD  pressure  tap  data  at  a  =  12° 
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CFD  pressure  tap  data  at  a 
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Figure  21-  CFD  pressure  tap  data  at  a  —  16°, 
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Figure  22.  CFD  pressure  tap  data  at  a  =  18° * 
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(a)  Lift  coefficient  (b)  Drag  coefficient 


(c)  Pitch-moment  coefficient 


Figure  23.  Unsteady  simulations  using  SARC-DDES  turbulence  model.  SARC-DDES  cases  ran  for  12,000 
iterations  and  coefficients  were  averaged  for  last  0,000  iterations. 
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